text:
"""
Generate the `src/services/{service}.jl` file.
"""
function _generate_high_level_wrapper(
    service_files::AbstractArray{ServiceFile}, auth::GitHub.OAuth2
)
    service_dir = joinpath(@__DIR__, "..", "services")

    # Remove old service files to ensure services that no longer exist are removed.
    for file in readdir(service_dir)
        path = joinpath(service_dir, file)
        if endswith(path, ".jl")
            rm(path)
        end
    end

    Threads.@threads for service_file in service_files
        service_name = service_file.name
        @info "Generating high level wrapper for $service_name"
        service = service_definition(service_file; auth=auth)

        service_name = lowercase(service["metadata"]["serviceId"])
        service_name = replace(service_name, ' ' => '_')
        operations = service["operations"]
        shapes = service["shapes"]

        protocol = service["metadata"]["protocol"]

        operations = sort!(
            _generate_high_level_definitions(service_name, protocol, operations, shapes)
        )

        service_path = joinpath(service_dir, "$service_name.jl")
        open(service_path, "w") do f
            println(
                f,
                """
             # This file is auto-generated by AWSMetadata.jl
             using AWS
             using AWS.AWSServices: $service_name
             using AWS.Compat
             using AWS.UUIDs
             """,
            )
            join(f, operations, "\n")
        end
    end
end

"""
Generate high-level definitions for the `service`.
All high-level definitions and documentation to be written into `services/{Service}.jl`
"""
function _generate_high_level_definitions(
    service_name::String, protocol::String, operations::AbstractDict, shapes::AbstractDict
)
    operation_definitions = String[]

    for (_, operation) in operations
        name = operation["name"]
        method = operation["http"]["method"]
        request_uri = operation["http"]["requestUri"]

        documentation = ""

        if haskey(operation, "documentation")
            documentation = _clean_documentation(operation["documentation"])
        end

        required_parameters = Dict{String,Any}()
        optional_parameters = Dict{String,Any}()

        if haskey(operation, "input")
            required_parameters, optional_parameters = _get_function_parameters(
                operation["input"]["shape"], shapes
            )
        end

        operation_definition = _generate_high_level_definition(
            service_name,
            protocol,
            name,
            method,
            request_uri,
            required_parameters,
            optional_parameters,
            documentation,
        )

        push!(operation_definitions, operation_definition)
    end

    return operation_definitions
end

"""
Generate the high-level definition for a services function.
"""
function _generate_high_level_definition(
    service_name::String,
    protocol::String,
    name::String,
    method::String,
    request_uri::String,
    required_parameters::AbstractDict,
    optional_parameters::AbstractDict,
    documentation::String,
)
    """
    Generate function definition for a service request given required, header and idempotent parameters.
    """
    function _generate_rest_operation_defintion(
        required_params::AbstractDict,
        optional_params::AbstractDict,
        function_name::String,
        service_name::String,
        method::String,
        request_uri::String,
    )
        request_uri = replace(request_uri, '{' => "\$(")  # Replace { with $(
        request_uri = replace(request_uri, '}' => ')')  # Replace } with )
        request_uri = replace(request_uri, '+' => "")  # Remove + from the request URI

        # Pre Julia-1.3 workaround
        req_keys = [replace(key, "-" => "_") for key in collect(keys(required_params))]

        required_params = filter(p -> (p[2]["location"] != "uri"), required_params)
        header_params = filter(p -> (p[2]["location"] == "header"), required_params)
        required_params = setdiff(required_params, header_params)

        idempotent_params = filter(p -> (p[2]["idempotent"]), optional_params)

        req_kv = ["\"$(p[1])\"=>$(replace(p[1], "-" => "_"))" for p in required_params]
        header_kv = ["\"$(p[1])\"=>$(replace(p[1], "-" => "_"))" for p in header_params]
        idempotent_kv = ["\"$(p[1])\"=>string(uuid4())" for p in idempotent_params]

        required_keys = !isempty(req_keys)
        headers = !isempty(header_params)
        idempotent = !isempty(idempotent_params)

        req_str = !isempty(req_kv) ? "Dict{String, Any}($(join(req_kv, ", "))" : ""
        params_str = if (!isempty(req_kv) || idempotent)
            "$(join(vcat(req_kv, idempotent_kv), ", "))"
        else
            ""
        end
        headers_str =
            headers ? "\"headers\"=>Dict{String, Any}($(join(header_kv, ", ")))" : ""
        params_headers_str = "Dict{String, Any}($(join([s for s in (params_str, headers_str) if !isempty(s)], ", ")))"

        formatted_function_name = _format_name(function_name)

        if required_keys && (idempotent || headers)
            return """
                $formatted_function_name($(join(req_keys, ", ")); aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", $params_headers_str; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                $formatted_function_name($(join(req_keys, ", ")), params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", Dict{String, Any}(mergewith(_merge, $params_headers_str, params)); aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                """
        elseif !required_keys && (idempotent || headers)
            return """
                $formatted_function_name(; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", $params_headers_str; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                $formatted_function_name(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", Dict{String, Any}(mergewith(_merge, $params_headers_str, params)); aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                """
        elseif required_keys && !isempty(req_kv)
            return """
                $formatted_function_name($(join(req_keys, ", ")); aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", $req_str); aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                $formatted_function_name($(join(req_keys, ", ")), params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", Dict{String, Any}(mergewith(_merge, $req_str), params)); aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                """
        elseif required_keys
            return """
                $formatted_function_name($(join(req_keys, ", ")); aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                $formatted_function_name($(join(req_keys, ", ")), params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                """
        else
            return """
                $formatted_function_name(; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                $formatted_function_name(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                """
        end
    end

    """
    Generate a JSON/Query high level definition.
    """
    function _generate_json_query_opeation_definition(
        required_params::AbstractDict,
        optional_params::AbstractDict,
        function_name::String,
        service_name::String,
    )
        req_keys = [replace(key, '-' => '_') for key in collect(keys(required_params))]

        idempotent_params = filter(p -> (p[2]["idempotent"]), optional_params)

        req_kv = ["\"$(p[1])\"=>$(replace(p[1], "-" => "_"))" for p in required_params]
        idempotent_kv = ["\"$(p[1])\"=>string(uuid4())" for p in idempotent_params]

        required = !isempty(req_kv)
        idempotent = !isempty(idempotent_kv)

        formatted_function_name = _format_name(function_name)

        if required && idempotent
            return """
                $formatted_function_name($(join(req_keys, ", ")); aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}($(join(req_kv, ", ")), $(join(idempotent_kv, ", "))); aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                $formatted_function_name($(join(req_keys, ", ")), params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}(mergewith(_merge, Dict{String, Any}($(join(req_kv, ", ")), $(join(idempotent_kv, ", "))), params)); aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                """
        elseif required
            return """
                $formatted_function_name($(join(req_keys, ", ")); aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}($(join(req_kv, ", "))); aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                $formatted_function_name($(join(req_keys, ", ")), params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}(mergewith(_merge, Dict{String, Any}($(join(req_kv, ", "))), params)); aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                """
        elseif idempotent
            return """
                $formatted_function_name(; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}($(join(idempotent_kv, ", "))); aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                $formatted_function_name(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}(mergewith(_merge, Dict{String, Any}($(join(idempotent_kv, ", "))), params)); aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                """
        else
            return """
                $formatted_function_name(; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$function_name\"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                $formatted_function_name(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = $service_name(\"$function_name\", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
                """
        end
    end

    """
    Generate the docstring for the `function_name`.
    """
    function _generate_docstring(
        function_name, documentation, required_parameters, optional_parameters
    )
        function_name = _format_name(function_name)
        args = join((_format_name(key) for (key, val) in required_parameters), ", ")
        maybejoin = isempty(args) ? "" : ", "
        operation_definition = """
            $(repeat('"', 3))
                $function_name($(args))
                $function_name($(args)$(maybejoin)params::Dict{String,<:Any})

            $(_wraplines(documentation))\n
            """

        # Add in the required parameters if applicable
        if !isempty(required_parameters)
            operation_definition *= "# Arguments\n"

            for (required_key, required_value) in required_parameters
                key = _format_name(required_key)
                operation_definition *= _wraplines(
                    "- `$key`: $(required_value["documentation"])"; delim="\n  "
                )
                operation_definition *= "\n"
            end

            operation_definition *= "\n"
        end

        # Add in the optional parameters if applicable
        if !isempty(optional_parameters)
            operation_definition *= """
                # Optional Parameters
                Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
                """

            for (optional_key, optional_value) in optional_parameters
                operation_definition *= _wraplines(
                    "- `\"$optional_key\"`: $(optional_value["documentation"])";
                    delim="\n  ",
                )
                operation_definition *= "\n"
            end
        end

        return operation_definition *= repeat('"', 3)
    end

    doc_string = _generate_docstring(
        name, documentation, required_parameters, optional_parameters
    )

    if protocol in ("json", "query", "ec2")
        function_string = _generate_json_query_opeation_definition(
            required_parameters, optional_parameters, name, service_name
        )
    elseif protocol in ("rest-json", "rest-xml")
        function_string = _generate_rest_operation_defintion(
            required_parameters,
            optional_parameters,
            name,
            service_name,
            method,
            request_uri,
        )
    else
        throw(
            ProtocolNotDefined(
                "$function_name is using a new protocol; $protocol which is not supported."
            ),
        )
    end

    return string(doc_string, '\n', function_string)
end


text:
# postprocessing.jl

using Plots

function plotSimulation(simulation, to::TimerOutput)
    # Plotting of 3D model
    @synctimeit to "plot result" begin
        grid = simulation.grid
        state = simulation.state

        ilayer = grid.nz

        nx = grid.nx
        ny = grid.ny
        nz = grid.nz

        x = grid.x
        y = grid.y
        z = grid.z

        # p1 = heatmap(x,y,state.h[:,:,ilayer]', fill = true, c=cgrad([:white,:green,:blue]))
        p1 = plot(x[0:nx + 1], y[0:ny + 1], state.hⁿ⁺¹[0:nx + 1,0:ny + 1,ilayer], seriestype=:wireframe, camera=(30, 60))
        # p2 = heatmap(state.u[:,nn,:]', fill = true, c=cgrad([:white,:green,:blue]))
        # p3 = heatmap(state.v[:,nn,:]', fill = true, c=cgrad([:white,:green,:blue]))
        # p4 = heatmap(state.w[:,nn,:]', fill = true, c=cgrad([:white,:green,:blue]))
        p2 = plot(x[0:nx + 1], state.hⁿ⁺¹[0:nx + 1,Int(ceil(grid.ny / 2)),ilayer])
        # p3 = plot(x[:],state.uⁿ⁺¹[1:nx,Int(grid.ny/2),ilayer])
        p3 = plot(z[0:nz + 1], state.wⁿ⁺¹[Int(ceil(grid.nx / 2)),Int(ceil(grid.nx / 2)),0:nz + 1])
        # Display the plots
        display(plot(p1, p2, p3, layout=(3, 1), size=(800, 1200)))
        # display(plot!(p1,p2,p3,p4, title=["water table h [m]" "horizontal velocity u [m/s]" "horizontal velocity v [m/s]" "vertical velocity w [m/s]"], layout=(4,1)))
        # savefig("plot.png")

    end # end timer
end

ext: jl


text:
atrisk_backedges = with_all_backedges(keys(becounter))

function atrisk_method(m::Method, atrisk_backedges)
    for mi in methodinstances(m)
        mi ∈ atrisk_backedges && return true
    end
    return false
end

function atrisk_triggers(m::Method, atrisk_instances)
    triggers = Set{Core.MethodInstance}()
    for mi in atrisk_instances
        if atrisk_method(m, all_backedges(mi))
            push!(triggers, mi)
        end
    end
    return triggers
end

# This removes MethodInstances that no one in their right mind should ever invalidate by specialization.
function remove_unlikely_methodinstances(list)
    out = Core.MethodInstance[]
    for mi in list
        mi = mi::Core.MethodInstance   # must have MethodInstance elements
        # All `continue` statements below omit the MethodInstance
        name = codename(mi.def.name)
        name ∈ (:invokelatest, :unwrap_unionall, :rewrap_unionall) && continue
        name ∈ (:print, :sprint) && length(mi.specTypes.parameters) - mi.def.nkw > 3 && continue
        # No one should ever specialize on notify or schedule's `val` argument
        name === :notify && !is_atrisk_type(mi.specTypes.parameters[2]) &&
            !any(is_atrisk_type, mi.specTypes.parameters[4:end]) && continue
        name === :schedule && !any(is_atrisk_type, mi.specTypes.parameters[2:end-1]) && continue
        # Add more removal-filters here

        # We've decided to keep it
        push!(out, mi)
    end
    return out
end

using Test

# Invalidating the code that loads packages leads to major slowdowns, especially if it happens repeatedly
# in a dependent chain of package loads. Ideally, we'd make this code-path "bulletproof".
for m in methods(Base.require)
    @test_broken isempty(remove_unlikely_methodinstances(atrisk_triggers(m, first.(miexp))))
    @test_broken isempty(remove_unlikely_methodinstances(atrisk_triggers(m, first.(mipriv))))
end

# Test overall number of atrisk MethodInstances and their average number of backedges
badexp = Set(remove_unlikely_methodinstances(first.(miexp)))
badcounts = filter(pr->pr.first ∈ badexp, miexp)
@test length(badcounts) < 1250  # original target 1000
if length(badcounts) < 1000 # 800
    @info "There are now only $(length(badcounts)) at-risk specializations of exported methods, consider dropping the threshold"
end
meancounts = sum(last.(badcounts))/length(badcounts)
@test meancounts < 32
if meancounts < 24
    @info "The mean number of at-risk backedges is now only $meancounts, consider dropping the threshold"
end

# Check for inference quality in specific functions.
# This is valid only for functions that should always return a particular type for any valid call of their methods.
function function_returns(@nospecialize(f), @nospecialize(typ); allow_missing_for_missing=true, minargs=0)
    for m in methods(f)
        sig = Base.unwrap_unionall(m.sig)
        for rt in Base.return_types(call_type(Base.unwrap_unionall(m.sig))...)
            rt <: typ && continue
            if allow_missing_for_missing && any(T->T===Missing, sig.parameters[2:end]) && rt === Missing
                continue
            end
            length(sig.parameters) - 1 < minargs && continue
            return false
        end
    end
    return true
end

# All the is* functions
# Not all of the broken cases have been checked carefully; it's possible some of these should return `Union{Bool,Missing}`
# or something.
@test_broken function_returns(isabspath, Bool)
@test function_returns(isabstracttype, Bool)
@test_broken function_returns(isapprox, Bool)
@test_broken function_returns(isascii, Bool)
# @test function_returns(isassigned, Bool)
@test function_returns(isbits, Bool)
@test function_returns(isbitstype, Bool)
@test function_returns(isblockdev, Bool)
@test function_returns(ischardev, Bool)
@test function_returns(iscntrl, Bool)
@test function_returns(isconcretetype, Bool)
@test function_returns(isconst, Bool)
@test function_returns(isdefined, Bool)
@test function_returns(isdigit, Bool)
@test function_returns(isdir, Bool)
@test function_returns(isdirpath, Bool)
@test_broken function_returns(isdisjoint, Bool)
@test function_returns(isdispatchtuple, Bool)
@test_broken function_returns(isempty, Bool)
@test_broken function_returns(isequal, Bool; minargs=2)
@test_broken function_returns(iseven, Bool)
@test function_returns(isexported, Bool)
@test function_returns(isfifo, Bool)
@test function_returns(isfile, Bool)
@test_broken function_returns(isfinite, Bool)
@test_broken function_returns(isinf, Bool)
@test_broken function_returns(isinteger, Bool)
@test function_returns(isinteractive, Bool)
@test_broken function_returns(isless, Bool)
@test function_returns(isletter, Bool)
@test function_returns(islink, Bool)
@test function_returns(islocked, Bool)
@test function_returns(islowercase, Bool)
@test_broken function_returns(ismarked, Bool)
@test function_returns(ismissing, Bool)
@test function_returns(ismount, Bool)
@test function_returns(ismutable, Bool)
@test function_returns(isnan, Bool)
@test function_returns(isnothing, Bool)
@test function_returns(isnumeric, Bool)
@test_broken function_returns(isodd, Bool)
@test_broken function_returns(isone, Bool)
@test_broken function_returns(isopen, Bool)
@test function_returns(ispath, Bool)
@test_broken function_returns(isperm, Bool)
@test_broken function_returns(ispow2, Bool)
@test function_returns(isprimitivetype, Bool)
@test function_returns(isprint, Bool)
@test function_returns(ispunct, Bool)
@test_broken function_returns(isreadable, Bool)
@test_broken function_returns(isreadonly, Bool)
@test_broken function_returns(isready, Bool)
@test_broken function_returns(isreal, Bool)
@test function_returns(issetequal, Bool)
@test function_returns(issetgid, Bool)
@test function_returns(issetuid, Bool)
@test function_returns(issocket, Bool)
@test_broken function_returns(issorted, Bool)
@test function_returns(isspace, Bool)
@test function_returns(issticky, Bool)
@test function_returns(isstructtype, Bool)
@test_broken function_returns(issubnormal, Bool)
@test_broken function_returns(issubset, Bool)
@test function_returns(istaskdone, Bool)
@test function_returns(istaskfailed, Bool)
@test function_returns(istaskstarted, Bool)
@test_broken function_returns(istextmime, Bool)
@test function_returns(isuppercase, Bool)
@test_broken function_returns(isvalid, Bool)
@test_broken function_returns(iswritable, Bool)
@test function_returns(isxdigit, Bool)
@test_broken function_returns(iszero, Bool)

@test function_returns(eof, Bool)

# Check that we never infer certain methodinstances
for f in (==, isequal, <, <=)
    for mi in methodinstances(f)
        if any(T->T<:Real, mi.specTypes.parameters)
            @test !is_atrisk_type(mi.specTypes)
        end
    end
end


text:
## bayes.jl
## (c) 2014, 2015 David A. van Leeuwen
##
## Attempt to implement a Bayesian approach to EM for GMMs, along the lines of
## Christopher Bishop's book, section 10.2.

## PLAN
## - convert stats collection to unnormalized stats and reduce results
## - optimize for speed

## Please note our pedantic use of the Greek symbols nu "ν" and alpha "α" (not to be confused with
## Latin v and a), and index "₀" and even superscript "⁻¹", which are both part of the identifiers
## e.g., W₀⁻¹ where others might write W0inv.

using SpecialFunctions: digamma

## initialize a prior with minimal knowledge
function GMMprior{T<:AbstractFloat}(d::Int, alpha::T, beta::T)
    m₀ = zeros(T, d)
    W₀ = eye(T, d)
    ν₀ = convert(T,d)
    GMMprior(alpha, beta, m₀, ν₀, W₀)
end
Base.copy(p::GMMprior) = GMMprior(p.α₀, p.β₀, copy(p.m₀), p.ν₀, copy(p.W₀))

## initialize from a GMM and nₓ, the number of points used to train the GMM.
function VGMM{T}(g::GMM{T}, π::GMMprior{T})
    nₓ = g.nx
    N = g.w * nₓ
    mx = g.μ
    if kind(g) == :diag
        S = covars(full(g))
    else
        S = covars(g)
    end
    α, β, m, ν, W = mstep(π, N, mx, S)
    hist = copy(g.hist)
    push!(hist, History("GMM converted to Variational GMM"))
    VGMM(g.n, g.d, π, α, β, m, ν, W, hist)
end
Base.copy(vg::VGMM) = VGMM(vg.n, vg.d, copy(vg.π), copy(vg.α), copy(vg.β),
                           copy(vg.m), copy(vg.ν), copy(vg.W), copy(vg.hist))

## W[k] really is chol(W_k, :U), use precision() to get it back
precision(c::AbstractTriangular) = c' * c
mylogdet(c::AbstractTriangular) = 2sum(log.(diag(c)))
mylogdet(m::Matrix) = logdet(m)

## sharpen VGMM to a GMM
## This currently breaks because my expected Λ are not positive definite
function GMM(vg::VGMM)
    w = vg.α / sum(vg.α)
    μ = vg.m
    Σ = similar(vg.W)
    for k=1:length(vg.W)
        Σ[k] = vg.W[k] * √vg.ν[k]
    end
    hist = copy(vg.hist)
    push!(hist, History("Variational GMM converted to GMM"))
    GMM(w, μ, Σ, hist, iround(sum(vg.α - vg.π.α₀)))
end

## m-step given prior and stats
function mstep{T}(π::GMMprior, N::Vector{T}, mx::Matrix{T}, S::Vector)
    ng = length(N)
    α = π.α₀ + N                # ng, 10.58
    ν = π.ν₀ + N + 1            # ng, 10.63
    β = π.β₀ + N                # ng, 10.60
    m = similar(mx)             # ng × d
    W = Array{eltype(FullCov{T})}(ng) # ng × (d*d)
    d = size(mx,2)
    limit = √ eps(eltype(N))
    W₀⁻¹ = inv(π.W₀)
    for k=1:ng
        if N[k] > limit
            m[k,:] = (π.β₀*π.m₀ + N[k]vec(mx[k,:])) ./ β[k] # 10.61 ## v0.5 arraymageddon
            Δ = vec(mx[k,:]) - π.m₀ ## v0.5 arraymageddon
            ## do some effort to keep the matrix positive definite
            third = π.β₀ * N[k] / (π.β₀ + N[k]) * (Δ * Δ') # guarantee symmety in Δ Δ'
            W[k] = chol(inv(cholfact(Symmetric(W₀⁻¹ + N[k]*S[k] + third)))) # 10.62
        else
            m[k,:] = zeros(d)
            W[k] = chol(eye(d))
        end
    end
    return α, β, m, ν, W
end

## this can be computed independently of the data
function expectations(vg::VGMM)
    ng, d = vg.n, vg.d
    Elogπ = digamma.(vg.α) .- digamma(sum(vg.α)) # 10.66, size ng
    ElogdetΛ = similar(Elogπ) # size ng
    for k in 1:ng
        ElogdetΛ[k] = sum(digamma.(0.5(vg.ν[k] .+ 1 - collect(1:d)))) + d*log(2) + mylogdet(vg.W[k]) # 10.65
    end
    return Elogπ, ElogdetΛ
end

## log(ρ_nk) from 10.46, start with a very slow implementation, and optimize it further
## This is as the heart of VGMM training, it should be super efficient, possibly
## at the cost of readability.
## Removing the inner loop over nₓ, replacing it by a matmul led to a major speed increase
## not leading to further speed increase for (Δ * vg.W[k]) .* Δ
## - c = Δ * chol(vg.W[k],:L), c .* c
## - Base.BLAS.symm('R', 'U', vg.ν[k], vg.W[k], Δ) .* Δ
function logρ(vg::VGMM, x::Matrix, ex::Tuple)
    (nₓ, d) = size(x)
    d == vg.d || error("dimension mismatch")
    ng = vg.n
    EμΛ = similar(x, nₓ, ng)    # nₓ × ng
    Δ = similar(x)              # nₓ × d
    for k in 1:ng
        ### d/vg.β[k] + vg.ν[k] * (x_i - m_k)' W_k (x_i = m_k) forall i
        ## Δ = (x_i - m_k)' W_k (x_i = m_k)
        xμTΛxμ!(Δ, x, vec(vg.m[k,:]), vg.W[k])
        EμΛ[:,k] = d/vg.β[k] .+ vg.ν[k] * sum(abs2, Δ, 2)
    end
    Elogπ, ElogdetΛ = ex
    (Elogπ + 0.5ElogdetΛ .- 0.5d*log(2π))' .- 0.5EμΛ
end

## r_nk 10.49, plus E[log p(Z|π)] 10.72 and E[log q(Z)] 10.75
function rnk(vg::VGMM, x::Matrix, ex::Tuple)
#    ρ = exp(logρ(g, x))
#    broadcast(/, ρ, sum(ρ, 2))
    logr = logρ(vg, x, ex)        # nₓ × ng
    broadcast!(-, logr, logr, logsumexp(logr, 2))
    ## we slowly replace logr by r, this is just cosmetic!
    r = logr
    ## ElogpZπ = sum(broadcast(*, r, Elogπ)) # E[log p(Z|π)] 10.72
    ElogpZπ = 0.                # E[log p(Z|π)] 10.72
    ElogqZ = 0.                 # E[log q(Z)] 10.75
    for k in 1:vg.n
        elπ = ex[1][k]          # Elogπ[k]
        for i = 1:size(logr,1)
            rr = exp(logr[i,k])
            ElogpZπ += rr * elπ
            ElogqZ += rr * logr[i,k]
            r[i,k] = rr         # cosmetic, r is identical to logr
        end
    end
    r, ElogpZπ + ElogqZ
end

## OK, also do stats.
## Like for the GMM, we return nₓ, (some value), zeroth, first, second order stats
## All return values can be accumulated, except r, which we need for
## lowerbound ElogpZπ and ElogqZ
function stats{T}(vg::VGMM, x::Matrix{T}, ex::Tuple)
    ng = vg.n
    (nₓ, d) = size(x)
    if nₓ == 0
        return 0, zero(RT), zeros(RT, ng), zeros(RT, ng, d), [zeros(RT, d,d) for k=1:ng]
    end
    r, ElogpZπqZ = rnk(vg, x, ex) # nₓ × ng
    N = vec(sum(r, 1))          # ng
    F = r' * x                  # ng × d
    ## S_k = sum_i r_ik x_i x_i'
    ## Sm = x' * hcat([broadcast(*, r[:,k], x) for k=1:ng]...)
    SS = similar(x, nₓ, d*ng)   # nₓ*nd*ng mem, nₓ*nd*ng multiplications
    for k in 1:ng
        for j=1:d for i=1:nₓ
            @inbounds SS[i,(k-1)*d+j] = r[i,k]*x[i,j]
        end end
    end
    Sm = x' * SS                # d * (d * ng) mem
    S = Matrix{T}[Sm[:,(k-1)*d+1:k*d] for k=1:ng]
    return nₓ, ElogpZπqZ, N, F, S
end

function stats(vg::VGMM, d::Data, ex::Tuple; parallel=false)
    if parallel
        r = dmap(x->stats(vg, x, ex), d)
        return reduce(+, r)
    else
        r = stats(vg, d[1], ex)
        for i=2:length(d)
            r += stats(vg, d[i], ex)
        end
        return r
    end
end

## trace(A*B) = sum(A' .* B)
function trAB{T1,T2}(A::Matrix{T1}, B::Matrix{T2})
    RT = promote_type(T1,T2)
    nr, nc = size(A)
    size(B) == (nc, nr) || error("Inconsistent matrix size")
    s=zero(RT)
    for i=1:nr for j=1:nc
        @inbounds s += A[i,j] * B[j,i]
    end end
    return s
end

## lower bound to the likelihood, using lots of intermediate results (10.2.2)
## ``We can straightforwardly evaluate the lower bound...''
function lowerbound(vg::VGMM, N::Vector, mx::Matrix, S::Vector,
                    Elogπ::Vector, ElogdetΛ::Vector, ElogpZπqZ)
    ## shorthands that make the formulas easier to read...
    ng, d = vg.n, vg.d
    α₀, β₀, ν₀, m₀, W₀  = vg.π.α₀, vg.π.β₀, vg.π.ν₀, vg.π.m₀, vg.π.W₀ # prior vars
    W₀⁻¹ = inv(W₀)
    α, β, m, ν, W = vg.α, vg.β, vg.m, vg.ν, vg.W # VGMM vars
    gaussians = 1:ng
    ## B.79
    logB(W, ν) = -0.5ν * (mylogdet(W) + d*log(2)) - d*(d-1)/4*log(π) - sum(lgamma.(0.5(ν + 1 - collect(1:d))))
    ## E[log p(x|Ζ,μ,Λ)] 10.71
    Elogll = 0.
    for k in gaussians
        Δ = vec(mx[k,:] - m[k,:])   # d ## v0.5 arraymageddon
        Wk = precision(W[k])
        Elogll += 0.5N[k] * (ElogdetΛ[k] - d/β[k] - d*log(2π)
                             - ν[k] * (trAB(S[k], Wk) + Wk * Δ ⋅ Δ)) # check chol efficiency
    end                        # 10.71
    ## E[log p(Z|π)] from rnk() 10.72
    Elogpπ = lgamma(ng*α₀) - ng*lgamma(α₀) - (α₀-1)sum(Elogπ) # E[log p(π)] 10.73
    ElogpμΛ = ng*logB(W₀, ν₀)   # E[log p(μ, Λ)] B.79
    for k in gaussians
        Δ = vec(m[k,:]) - m₀        # d ## v0.5 arraymageddon
        Wk = precision(W[k])
        ElogpμΛ += 0.5(d*log(β₀/(2π)) + (ν₀-d)ElogdetΛ[k] - d*β₀/β[k]
                       -β₀*ν[k] * Wk * Δ ⋅ Δ - ν[k]*trAB(W₀⁻¹, Wk))
    end                         # 10.74
    ## E[log q(Z)] from rnk() 10.75, combined with E[log p(Z|π)]
    Elogqπ = sum((α.-1).*Elogπ) + lgamma(sum(α)) - sum(lgamma.(α)) # E[log q(π)] 10.76
    ## E[log q(μ,Λ)] 10.77
    ElogqμΛ = 0.
    for k in gaussians
        H = -logB(W[k],ν[k]) - 0.5(ν[k]-d-1)ElogdetΛ[k] + ν[k]*d/2 # H[q(Λ)] B.82
        ElogqμΛ += 0.5(ElogdetΛ[k] + d*log(β[k]/(2π)) - d) - H
    end                         # 10.77
    return Elogll + ElogpZπqZ + Elogpπ + ElogpμΛ + Elogqπ + ElogqμΛ
end

## note: for matrix argument, Gaussian index must run down!
rmdisfunct(v::Vector, keep) = v[keep]
rmdisfunct(m::Matrix, keep) = m[keep, :]

## do exactly one update step for the VGMM, and return an estimate for the lower bound
## of the log marginal probability p(x)
function emstep!(vg::VGMM, x::DataOrMatrix)
    ## E-like step
    Elogπ, ElogdetΛ = expectations(vg)
    ## N, mx, S, r = threestats(vg, x, (Elogπ, ElogdetΛ))
    nₓ, ElogZπqZ, N, F, S = stats(vg, x, (Elogπ, ElogdetΛ))
    mx = F ./ N
    for k = 1:vg.n
        mk = vec(mx[k,:]) ## v0.5 arraymageddon
        S[k] = S[k] / N[k] - mk * mk'
    end
    ## remove defunct Gaussians
    keep = N .> √ eps(eltype(N))
    n = sum(keep)
    if n < vg.n
        N, mx, S, Elogπ, ElogdetΛ = map(x->rmdisfunct(x,keep), (N, mx, S, Elogπ, ElogdetΛ))
        vg.α, vg.β, vg.m, vg.ν, vg.W = map(x->rmdisfunct(x,keep), (vg.α, vg.β, vg.m, vg.ν, vg.W))
        vg.n = n
        addhist!(vg, @sprintf("dropping number of Gaussions to %d",n))
    end
    ## then compute the lowerbound
    L = lowerbound(vg, N, mx, S, Elogπ, ElogdetΛ, ElogZπqZ)
    ## and finally compute M-like step
    vg.α, vg.β, vg.m, vg.ν, vg.W = mstep(vg.π, N, mx, S)
    L, nₓ
end

## This is called em!, but it is not really expectation maximization I think
function em!(vg::VGMM, x; nIter=50)
    L = Float64[]
    nₓ = 0
    for i in 1:nIter
        lb, nₓ = emstep!(vg, x)
        push!(L, lb)
        if i > 1 && isapprox(L[i], L[i-1], rtol=0)
            nIter = i
            break
        end
        addhist!(vg, @sprintf("iteration %d, lowerbound %f", i, last(L)/nₓ/vg.d))
    end
    L ./= nₓ*vg.d
    addhist!(vg, @sprintf("%d variational Bayes EM-like iterations using %d data points, final lowerbound %f", nIter, nₓ, last(L)))
    L
end

## Not used

if false
## Wishart distribution {\cal W}(Λ, W, ν).
function Wishart(Λ::Matrix, W::Matrix, ν::Float64)
    ## check
    issym(W) || error("W must be symmetric")
    d = size(W,1)
    (d,d) == size(Λ) || error("Λ must be same size as W")
    ## norm
    B = det(W)^(-0.5ν) / (2^(0.5ν*d) * π^(0.25d*(d-1)))
    for i=1:d
        B /= gamma(0.5(ν+1-i))
    end
#    invΛ = inv(Λ)
    ex = -0.5dot(inv(W),Λ)
    B * det(Λ)^(0.5(ν-d-1)) * exp(ex)
end

function Gaussian(x::Vector, μ::Vector, Σ::Matrix)
    d = length(μ)
    length(x) == d || error("Wrong dimension x")
    size(Σ) = (d,d) || error("Inconsistent size Σ")
    norm = inv((2π)^0.5d * sqrt(det(Σ)))
    ex = -0.5(x-μ)' * inv(Σ) * (x-μ)
    norm * exp(ex)
end

function GaussianWishart(μ::Vector, Λ::Matrix, μ₀::Vector, β::Float64, W::Matrix, ν::Float64)
    Gaussian(μ, μ₀, inv(β*Λ)) * Wishart(Λ, W, ν)
end

end


text:
module TestUtils

import Random
import JuMP
using Test: @test

export test_inverse_solution, test_data_fidelity, add_noise

function noisy_observation(observation_model, x; rng = Random.MersenneTwister(1))
    ŷ = observation_model.expected_observation(x)
    ŷ + randn(rng, size(ŷ)) .* observation_model.σ
end

function test_inverse_solution(weights_est, weights_true; atol = 1e-2, verbose = false)
    w_total_est = sum(weights_est)
    w_total_true = sum(weights_true)

    for k in keys(weights_true)
        verbose && @info k
        @test isapprox(weights_est[k] / w_total_est, weights_true[k] / w_total_true; atol = atol)
    end
end

function test_data_fidelity(inverse_model, observation_model, x̂, y)
    ŷ = observation_model.expected_observation(x̂)
    ê_sq = sum((ŷ .- y) .^ 2)
    @test JuMP.objective_value(inverse_model) <= ê_sq + 1e-2
end

end


text:
# src path
src_path = "C:\\Users\\hidee\\Desktop\\git\\general_2d_NS_Chemical-N2-N-\\src_c\\"

# main (変更しないこと)
src_read="read_grid.jl"
include(src_path*src_read)
src_read="read_para.jl"
include(src_path*src_read)
src_read="output.jl"
include(src_path*src_read)

src_read="advection_term.jl"
include(src_path*src_read)
src_read="boundary.jl"
include(src_path*src_read)
src_read="converge.jl"
include(src_path*src_read)
src_read="misc.jl"
include(src_path*src_read)
src_read="rhs.jl"
include(src_path*src_read)
src_read="value_setup.jl"
include(src_path*src_read)
src_read="viscos.jl"
include(src_path*src_read)
src_read="chemical_constant.jl"
include(src_path*src_read)
src_read="chemical_term.jl"
include(src_path*src_read)
src_read="setup.jl"
include(src_path*src_read)

src_read="cal_time_step.jl"
include(src_path*src_read)
src_main="main.jl"
include(src_path*src_main)


text:
f = open(join([version_dir, "optim_linesearch_benchmark.csv"], "/"), "w")
write(f, join(["Problem", "Optimizer", "Converged", "Time", "Minimum", "Iterations", "f_calls", "g_calls", "h_call", "f_hat", "f_error", "x_error"], ","))
write(f, "\n")

linesearch_solver_names = ["Accelerated Gradient Descent",
                           "BFGS",
                           "L-BFGS",
                           "Conjugate Gradient",
                           "Momentum Gradient Descent",
                           "Gradient Descent",
                           "Newton"]
linesearch_solver_names = ["Newton",]

linesearch_solvers =[AcceleratedGradientDescent,
                            BFGS,
                            LBFGS,
                            ConjugateGradient,
                            MomentumGradientDescent,
                            GradientDescent,
                            Newton]
linesearch_solvers = [Newton, ]

default_linesearches = [(LineSearches.BackTracking(order = 2), "BackTracking (quadratic)"),
                        (LineSearches.BackTracking(order = 3), "BackTracking (qubic)"),
                        (LineSearches.HagerZhang(), "HagerZhang"),
                        (LineSearches.MoreThuente(), "MoreThuente"),
                        (LineSearches.Static(alpha=0.1), "Static (alpha = 0.1)"),
                        (LineSearches.Static(), "Static (alpha = 1)"),
                        (LineSearches.StrongWolfe(), "StrongWolfe")]

@showprogress 1 "Benchmarking..." for (name, problem) in Optim.UnconstrainedProblems.examples
    for (i,  algorithm) in enumerate(linesearch_solvers)
        for ls in default_linesearches
            print_debuginfo && @show name, problem , algorithm
            try
            # Force compilation and obtain results
            results = optimize(problem.f, problem.g!, problem.h!,
                               problem.initial_x, algorithm(linesearch=ls[1]), Optim.Options(g_tol = 1e-16))
            # Run each algorithm n times
            n = 10
            if algorithm == ParticleSwarm() && name == "Large Polynomial"
                n = 1
            end
            # Estimate run time in seconds
            run_time = minimum([@elapsed optimize(problem.f, problem.g!, problem.h!,
                                   problem.initial_x,
                                   algorithm(linesearch=ls[1]), Optim.Options(g_tol = 1e-16)) for nn = 1:n])

            # Count iterations
            iterations = results.iterations

            # Print out results.
            write(f, join([problem.name,
                           linesearch_solver_names[i]*" with "*ls[2],
                           Optim.converged(results),
                           run_time,
                           Optim.minimum(results),
                           iterations,
                           Optim.f_calls(results),
                           Optim.g_calls(results),
                           Optim.h_calls(results),
                           problem.f(problem.solutions),
                           Optim.minimum(results)-problem.f(problem.solutions),
                           norm(Optim.minimizer(results)-problem.solutions, Inf)], ","))
            write(f, "\n")
            catch m
                println(m)
                write(f, join([problem.name,
                               linesearch_solver_names[i]*" with "*ls[2],
                               "false",
                               "Inf",
                               "Inf",
                               "Inf",
                               "Inf",
                               "Inf",
                               "Inf",
                               "Inf",
                               "Inf",
                               "Inf"], ","))
                write(f, "\n")
            end
        end
    end
end
close(f)


text:
using HorizonSideRobots

function corner(r::Robot)
    n=0
    w=0
    while !isborder(r,Nord)
        move!(r,Nord)
        n+=1
    end
    while !isborder(r,West)
        move!(r,West)
        w+=1
    end
    all(r)
    while !isborder(r,Nord)
        move!(r,Nord)
    end
    while !isborder(r,West)
        move!(r,West)
    end
    for a in 1:w
        move!(r,Ost)
    end
    for a in 1:n
        move!(r,Sud)
    end
end

function all(r::Robot)
    while !isborder(r,Sud)
        for side in [Ost,West]
            while !isborder(r,side)
                putmarker!(r)
                move!(r,side)
            end
            putmarker!(r)
            move!(r,Sud)
        end
    end
    while !isborder(r,Ost)
        putmarker!(r)
        move!(r,Ost)
    end
    putmarker!(r)
end

corner(r)

max_stars_repo_name: Cheburkova/mirea-progs

max_issues_repo_name: Cheburkova/mirea-progs

max_forks_repo_name: Cheburkova/mirea-progs


text:
using CUDAapi

using Test

@testset "library types" begin
    @test CUDAapi.PATCH_LEVEL == CUDAapi.libraryPropertyType(2)
    @test CUDAapi.C_32U == CUDAapi.cudaDataType(Complex{UInt32})
end

@testset "properties" begin
    @test !CUDAapi.gcc_supported(v"5.0", v"5.5")
    @test CUDAapi.gcc_supported(v"5.0", v"8.0")
    CUDAapi.devices_for_cuda(v"8.0")
    CUDAapi.devices_for_llvm(v"5.0")
    CUDAapi.isas_for_cuda(v"8.0")
    CUDAapi.isas_for_llvm(v"5.0")
end

# helper macro to test for non-nothingness
macro test_something(ex)
    quote
        local rv = $(esc(ex))
        @test rv !== nothing
        rv
    end
end

@testset "discovery" begin
    find_binary([Sys.iswindows() ? "CHKDSK" : "true"])
    find_library([Sys.iswindows() ? "NTDLL" : "c"])

    dirs = find_toolkit()
    @test !isempty(dirs)

    ver = find_toolkit_version(dirs)

    @testset "CUDA tools and libraries" begin
        @test_something find_cuda_binary("nvcc", dirs)
        @test_something find_cuda_library("cudart", dirs)
        @test_something find_cuda_library("nvtx", dirs)
        @test_something find_libdevice([v"3.0"], dirs)
        @test_something find_libcudadevrt(dirs)
        @test_something find_toolchain(dirs)
        @test_something find_toolchain(dirs, ver)
    end

    @testset "host tools and libraries" begin
        @test_something find_host_compiler()
        @test_something find_host_compiler(ver)
    end

    if haskey(ENV, "CI")
        # CI deals with plenty of CUDA versions, which makes discovery tricky.
        # dump a relevant tree of files to help debugging
        function traverse(dir, level=0)
            for entry in readdir(dir)
                print("  "^level)
                path = joinpath(dir, entry)
                if isdir(path)
                    println("└ $entry")
                    traverse(path, level+1)
                else
                    println("├ $entry")
                end
            end
        end
        for dir in dirs
            println("File tree of toolkit directory $dir:")
            traverse(dir)
        end
    end
end

@testset "availability" begin
    @test isa(has_cuda(), Bool)
    @test isa(has_cuda_gpu(), Bool)
end


text:
export PiecewiseConstants

"A dictionary of piecewise constant functions on `[0,1]`."
struct PiecewiseConstants{T} <: Dictionary1d{T,T}
    n       ::  Int
end

PiecewiseConstants(n = 10) = PiecewiseConstants{Float64}(n)

## Basic interface

# size should return a tuple (length is defined automatically in terms of size)
size(h::PiecewiseConstants) = (h.n,)

# By implementing `similar`, we support promotion of the domain type, as well
# as resizing of the dictionary.
# If resizing is not possible, one can add a statement such as
# "@assert n == length(dict)", which will cause an error upon any user attempt to resize.
# Same for promotion.
similar(::PiecewiseConstants, ::Type{S}, n::Int) where {S} = PiecewiseConstants{S}(n)

# Evaluation of the basis function:
unsafe_eval_element(h::PiecewiseConstants{T}, idx::Int, x) where {T} = one(T)

# The interpolation_grid returns a set of points such that the interpolation
# problem is uniquely solvable. This is useful for function approximation.
interpolation_grid(h::PiecewiseConstants) = MidpointEquispacedGrid(length(h), support(h))

# The support of our dictionary is the half-open interval [0,1)
support(h::PiecewiseConstants{T}) where {T} = Interval{:closed,:open,T}(0,1)

# We also implement the support of each basis function. By default, this would
# be the default of the dictionary itself.
support(h::PiecewiseConstants{T}, idx::Int) where {T} =
    Interval{:closed,:open,T}(0,one(T)/h.n) .+ (one(T)/h.n*(idx-1))

## Orthogonality properties

# The measure of a dictionary is used when computing projections (inner products
# of a function with a basis function)
measure(h::PiecewiseConstants{T}) where {T} = FourierWeight{T}()

isorthogonal(h::PiecewiseConstants, measure::FourierWeight) = true

# Is our dictionary orthogonal with respect to a discrete measure? We check all points.
function isorthogonal(h::PiecewiseConstants, μ::DiscreteWeight)
    for (i,xi) in zip(ordering(h),points(μ))
        if xi ∉ support(h, i)
            return false
        end
    end
    return true
end


text:
# Pre-defined C-function API of the dynamically loaded library of openBuildNet

libobnext = C_NULL

function load_obnlibext()
    global libobnext
    if libobnext == C_NULL
        libobnext = Libdl.dlopen("libobnext-mqtt", Libdl.RTLD_GLOBAL | Libdl.RTLD_LAZY)

        if libobnext == C_NULL
            error("Failed to load the openBuildNet External Interface Library (libobnext-mqtt).")
        end
    end
end

load_obnlibext()

# The C API exported from the library
# See the header file obnnode_ext.h for details.

# Mask type to specify which blocks are to be triggered
typealias OBNUpdateMask UInt64

# Simulation time type
typealias OBNSimTimeType Int64

# Event type (C enum) and possible values
typealias OBNEI_EventType Cuint
const OBNEI_Event_INIT = convert(OBNEI_EventType, 0)
const OBNEI_Event_Y = convert(OBNEI_EventType, 1)
const OBNEI_Event_X = convert(OBNEI_EventType, 2)
const OBNEI_Event_TERM = convert(OBNEI_EventType, 3)
const OBNEI_Event_RCV = convert(OBNEI_EventType, 4)
const OBNEI_Event_RESTART = convert(OBNEI_EventType, 5)

# Type to pass arguments of an event (C struct)
immutable OBNEI_EventArg
  mask::OBNUpdateMask
  index::Csize_t
end

# Node interface
const _api_createOBNNode = Libdl.dlsym(libobnext::Ptr{Void}, :createOBNNode)
const _api_deleteOBNNode = Libdl.dlsym(libobnext::Ptr{Void}, :deleteOBNNode)
const _api_nodeStopSimulation = Libdl.dlsym(libobnext::Ptr{Void}, :nodeStopSimulation)
const _api_nodeRequestStopSimulation = Libdl.dlsym(libobnext::Ptr{Void}, :nodeRequestStopSimulation)

const _api_nodeIsStopped = Libdl.dlsym(libobnext::Ptr{Void}, :nodeIsStopped)
const _api_nodeIsError = Libdl.dlsym(libobnext::Ptr{Void}, :nodeIsError)
const _api_nodeIsRunning = Libdl.dlsym(libobnext::Ptr{Void}, :nodeIsRunning)

const _api_nodeSimulationTime = Libdl.dlsym(libobnext::Ptr{Void}, :nodeSimulationTime)
const _api_nodeSimulationTimeTicks = Libdl.dlsym(libobnext::Ptr{Void}, :nodeSimulationTimeTicks)
const _api_nodeWallClockTime = Libdl.dlsym(libobnext::Ptr{Void}, :nodeWallClockTime)
const _api_nodeTimeUnit = Libdl.dlsym(libobnext::Ptr{Void}, :nodeTimeUnit)

const _api_simGetPortEvent = Libdl.dlsym(libobnext::Ptr{Void}, :simGetPortEvent)
const _api_simRunStep = Libdl.dlsym(libobnext::Ptr{Void}, :simRunStep)
const _api_simSetEventResult = Libdl.dlsym(libobnext::Ptr{Void}, :simSetEventResult)
const _api_simRequestFutureUpdate = Libdl.dlsym(libobnext::Ptr{Void}, :simRequestFutureUpdate)

# Port interface
const _api_createInputPort = Libdl.dlsym(libobnext::Ptr{Void}, :createInputPort)
const _api_createOutputPort = Libdl.dlsym(libobnext::Ptr{Void}, :createOutputPort)
const _api_inputPending = Libdl.dlsym(libobnext::Ptr{Void}, :inputPending)
const _api_portConnect = Libdl.dlsym(libobnext::Ptr{Void}, :portConnect)
const _api_outputSendSync = Libdl.dlsym(libobnext::Ptr{Void}, :outputSendSync)
const _api_portInfo = Libdl.dlsym(libobnext::Ptr{Void}, :portInfo)
const _api_portEnableRcvEvent = Libdl.dlsym(libobnext::Ptr{Void}, :portEnableRcvEvent)

# Scalar inputs
const _api_inputScalarDoubleGet = Libdl.dlsym(libobnext::Ptr{Void}, :inputScalarDoubleGet)
const _api_inputScalarBoolGet = Libdl.dlsym(libobnext::Ptr{Void}, :inputScalarBoolGet)
const _api_inputScalarInt32Get = Libdl.dlsym(libobnext::Ptr{Void}, :inputScalarInt32Get)
const _api_inputScalarInt64Get = Libdl.dlsym(libobnext::Ptr{Void}, :inputScalarInt64Get)
const _api_inputScalarUInt32Get = Libdl.dlsym(libobnext::Ptr{Void}, :inputScalarUInt32Get)
const _api_inputScalarUInt64Get = Libdl.dlsym(libobnext::Ptr{Void}, :inputScalarUInt64Get)

# Vector inputs
const _api_inputVectorDoubleGet = Libdl.dlsym(libobnext::Ptr{Void}, :inputVectorDoubleGet)
const _api_inputVectorDoubleRelease = Libdl.dlsym(libobnext::Ptr{Void}, :inputVectorDoubleRelease)

const _api_inputVectorBoolGet = Libdl.dlsym(libobnext::Ptr{Void}, :inputVectorBoolGet)
const _api_inputVectorBoolRelease = Libdl.dlsym(libobnext::Ptr{Void}, :inputVectorBoolRelease)

const _api_inputVectorInt32Get = Libdl.dlsym(libobnext::Ptr{Void}, :inputVectorInt32Get)
const _api_inputVectorInt32Release = Libdl.dlsym(libobnext::Ptr{Void}, :inputVectorInt32Release)

const _api_inputVectorInt64Get = Libdl.dlsym(libobnext::Ptr{Void}, :inputVectorInt64Get)
const _api_inputVectorInt64Release = Libdl.dlsym(libobnext::Ptr{Void}, :inputVectorInt64Release)

const _api_inputVectorUInt32Get = Libdl.dlsym(libobnext::Ptr{Void}, :inputVectorUInt32Get)
const _api_inputVectorUInt32Release = Libdl.dlsym(libobnext::Ptr{Void}, :inputVectorUInt32Release)

const _api_inputVectorUInt64Get = Libdl.dlsym(libobnext::Ptr{Void}, :inputVectorUInt64Get)
const _api_inputVectorUInt64Release = Libdl.dlsym(libobnext::Ptr{Void}, :inputVectorUInt64Release)

# Matrix inputs
const _api_inputMatrixDoubleGet = Libdl.dlsym(libobnext::Ptr{Void}, :inputMatrixDoubleGet)
const _api_inputMatrixDoubleRelease = Libdl.dlsym(libobnext::Ptr{Void}, :inputMatrixDoubleRelease)

const _api_inputMatrixBoolGet = Libdl.dlsym(libobnext::Ptr{Void}, :inputMatrixBoolGet)
const _api_inputMatrixBoolRelease = Libdl.dlsym(libobnext::Ptr{Void}, :inputMatrixBoolRelease)

const _api_inputMatrixInt32Get = Libdl.dlsym(libobnext::Ptr{Void}, :inputMatrixInt32Get)
const _api_inputMatrixInt32Release = Libdl.dlsym(libobnext::Ptr{Void}, :inputMatrixInt32Release)

const _api_inputMatrixInt64Get = Libdl.dlsym(libobnext::Ptr{Void}, :inputMatrixInt64Get)
const _api_inputMatrixInt64Release = Libdl.dlsym(libobnext::Ptr{Void}, :inputMatrixInt64Release)

const _api_inputMatrixUInt32Get = Libdl.dlsym(libobnext::Ptr{Void}, :inputMatrixUInt32Get)
const _api_inputMatrixUInt32Release = Libdl.dlsym(libobnext::Ptr{Void}, :inputMatrixUInt32Release)

const _api_inputMatrixUInt64Get = Libdl.dlsym(libobnext::Ptr{Void}, :inputMatrixUInt64Get)
const _api_inputMatrixUInt64Release = Libdl.dlsym(libobnext::Ptr{Void}, :inputMatrixUInt64Release)

# Scalar outputs
const _api_outputScalarDoubleSet = Libdl.dlsym(libobnext::Ptr{Void}, :outputScalarDoubleSet)
const _api_outputScalarBoolSet = Libdl.dlsym(libobnext::Ptr{Void}, :outputScalarBoolSet)
const _api_outputScalarInt32Set = Libdl.dlsym(libobnext::Ptr{Void}, :outputScalarInt32Set)
const _api_outputScalarInt64Set = Libdl.dlsym(libobnext::Ptr{Void}, :outputScalarInt64Set)
const _api_outputScalarUInt32Set = Libdl.dlsym(libobnext::Ptr{Void}, :outputScalarUInt32Set)
const _api_outputScalarUInt64Set = Libdl.dlsym(libobnext::Ptr{Void}, :outputScalarUInt64Set)

# Vector outputs
const _api_outputVectorDoubleSet = Libdl.dlsym(libobnext::Ptr{Void}, :outputVectorDoubleSet)
const _api_outputVectorBoolSet = Libdl.dlsym(libobnext::Ptr{Void}, :outputVectorBoolSet)
const _api_outputVectorInt32Set = Libdl.dlsym(libobnext::Ptr{Void}, :outputVectorInt32Set)
const _api_outputVectorInt64Set = Libdl.dlsym(libobnext::Ptr{Void}, :outputVectorInt64Set)
const _api_outputVectorUInt32Set = Libdl.dlsym(libobnext::Ptr{Void}, :outputVectorUInt32Set)
const _api_outputVectorUInt64Set = Libdl.dlsym(libobnext::Ptr{Void}, :outputVectorUInt64Set)

# Matrix outputs
const _api_outputMatrixDoubleSet = Libdl.dlsym(libobnext::Ptr{Void}, :outputMatrixDoubleSet)
const _api_outputMatrixBoolSet = Libdl.dlsym(libobnext::Ptr{Void}, :outputMatrixBoolSet)
const _api_outputMatrixInt32Set = Libdl.dlsym(libobnext::Ptr{Void}, :outputMatrixInt32Set)
const _api_outputMatrixInt64Set = Libdl.dlsym(libobnext::Ptr{Void}, :outputMatrixInt64Set)
const _api_outputMatrixUInt32Set = Libdl.dlsym(libobnext::Ptr{Void}, :outputMatrixUInt32Set)
const _api_outputMatrixUInt64Set = Libdl.dlsym(libobnext::Ptr{Void}, :outputMatrixUInt64Set)

# Binary ports
const _api_inputBinaryGet = Libdl.dlsym(libobnext::Ptr{Void}, :inputBinaryGet)
const _api_inputBinaryRelease = Libdl.dlsym(libobnext::Ptr{Void}, :inputBinaryRelease)
const _api_outputBinarySet = Libdl.dlsym(libobnext::Ptr{Void}, :outputBinarySet)

# Misc
const _api_lastErrorMessage = Libdl.dlsym(libobnext::Ptr{Void}, :lastErrorMessage)
function lastErrorMessage()
  val = ccall(_api_lastErrorMessage, Ptr{UInt8}, ())
  if val == C_NULL
    error("Could not get OBN's last error message.")
  end
  bytestring(val)
end

const _api_maxUpdateID = Libdl.dlsym(libobnext::Ptr{Void}, :maxUpdateID)
max_blockid() = ccall(_api_maxUpdateID, Cint, ())

# Convert a list of block IDs into an update mask
function update_mask(blks...)
  mask = OBNUpdateMask(0)
  maxid = max_blockid()

  for id in blks
    iid = UInt64(id)  # Convert to unsigned int
    @assert iid<=maxid "Invalid block ID."

    mask |= (1 << iid)
  end
  mask
end


text:
function run_cross_validation(pm::PredictionModel, inference_config::InferenceConfiguration, model_config::Union{Nothing, ModelConfiguration}; num_folds::Int64 = 10)
    if !inference_config.do_cross_validation
        throw(ArgumentError("run_cross_validation called without proper inference configuration, make sure the do_cross_validation flag is set to true"))
    end
    # running cross validation means, we are discarding the original training set/test set designations
    num_cl = length(pm.data.cell_lines)
    cls_per_fold = Int64(floor(num_cl/num_folds))

    for i in 1:num_folds
        log_progress(pm.data.internal_id, "starting computation of fold $i in CV")
        log_message("starting fold $i")
        start_idx = (i-1)*cls_per_fold + 1
        end_idx = i == num_folds ? num_cl : i * cls_per_fold

        # mark cell lines as test and training set
        for (cl_idx, cl) in enumerate(collect(values(pm.data.cell_lines)))
            if start_idx ≤ cl_idx ≤ end_idx
                cl.in_test_set = true
            else
                cl.in_test_set = false
            end
        end
        inference_config.fold_num = i
        log_message("reshuffled training/test set and set fold to $i")
        run_model(pm, inference_config, model_config)
    end
end

function run_model(pm::PredictionModel, inference_config::InferenceConfiguration, model_config::Union{Nothing, ModelConfiguration})
    # things to do only once (or once per fold)
    # for BMTMKL, this computes kernels and cross kernels depending on which cell lines
    # are in the training and test set
    log_message("running post_init!")
    post_init!(pm)
    if model_config != nothing update_model_config!(pm, model_config) end
    # result dir is the target directory amended by the fold
    result_dir = inference_config.do_cross_validation ? joinpath(inference_config.target_dir, string(inference_config.fold_num)) : inference_config.target_dir
    mkpath(result_dir)

    if inference_config.compute_wpc_index
        # preparation for wpc index computations

        # 1) create files with test responses and training responses respectively (for this fold)
        # training
        training_response_file = create_response_file(pm, false, result_dir, "training_response.txt")
        # test
        test_response_file = create_response_file(pm, true, result_dir, "test_response.txt")

        # 2) create random predictions from these
        script_name = joinpath(PROJECT_ROOT, "analytics", "create_random_predictions.pl")
        rand_pred_dir = joinpath(result_dir, "random_predictions")
        mkpath(rand_pred_dir)
        create_rand_pred_cmd = `perl $script_name $training_response_file $test_response_file $rand_pred_dir 100`
        run(create_rand_pred_cmd)

        # 3) calculate z-scores
        script_name = joinpath(PROJECT_ROOT, "analytics", "calculate_zscore_weights.pl")
        sd_file = joinpath(PROJECT_ROOT, "analytics", "pooled_sd.txt")
        gold_standard_file = joinpath(PROJECT_ROOT, "analytics", "gold_standard.csv")
        zscore_file = joinpath(result_dir, "zscores.txt")
        calc_zscore_cmd = `perl $script_name $rand_pred_dir $sd_file $test_response_file $gold_standard_file $zscore_file`
        run(calc_zscore_cmd)

        # 4) with random prediction and zscores, compute the wpc index
        script_name = joinpath(PROJECT_ROOT, "analytics", "weighted_average_concordance_index.pl")
        mkpath(joinpath(result_dir, "wpc"))
        wpc_target_file = joinpath(result_dir, "wpc", "gold_standard_wpc.txt")

        wpc_summary_file = joinpath(inference_config.target_dir, "wpc_summary.txt")

        file = open(wpc_summary_file, "w")
        @printf(file, "ranking file\taverage probabilistic c-index\tweighted average probabilistic c-index\twpc p-value\n")
        close(file)

        wpc_script_gs_cmd = `perl $script_name $sd_file $test_response_file $zscore_file $gold_standard_file $wpc_target_file $wpc_summary_file`
        run(wpc_script_gs_cmd)
    end
    all_model_configs = ModelConfiguration[]
    # if doing gridsearch, create model configs to use
    if inference_config.do_gridsearch
        gamma_dist_alphas = [250., 500., 750., 1000.]
        alpha_beta_ratios = [.1, .5, 1., 2.5, 5.]
        normal_vars = [.1, .5, 1., 2.]
        # normal_means = [1.]
        for alpha in gamma_dist_alphas, ratio in alpha_beta_ratios, v in normal_vars
            mc = ModelConfiguration(pm, alpha, alpha/ratio, 1., v)
            push!(all_model_configs, mc)
        end
        log_message("generated gridsearch configs ($(length(all_model_configs)))")
    #  else just use the model config as provided in the call
    else
        model_config == nothing ?
            throw(ArgumentError("not doing a grid search but no model configuration specified")) :
            push!(all_model_configs, model_config)
    end
    all_errors = Vector{String}(undef, length(all_model_configs))

    # do the actual inference sweep
    cnt = Base.Threads.Atomic{Int64}(0)
    Threads.@threads for i in 1:length(all_model_configs)
        mc = all_model_configs[i]
        # TODO: this assumes we are using the same values for all params, maybe change later
        # same for the generic filenam below
        alpha = mc.parameters["α_ɣ"]
        beta = mc.parameters["β_ɣ"]
        mu = mc.parameters["μ_e"]
        v = mc.parameters["σ_e"]
        Base.Threads.atomic_add!(cnt, 1)
        log_message("parameter setting $(cnt[]) (of $(length(all_model_configs))): alpha=$alpha, beta=$beta mu=$mu var=$v")
        # TODO: run each model multiple times (e.g. 10) and collect prediction results
        # report mean and variance of those predictions
        # log_message("calling parameter_inference method")
        (lls, errs, test_errs, params, convergence) = inference(pm, mc, inference_config = inference_config)
        log_message("inference stats: likelihood=$(lls[end]) train_error=$(errs[end]) test_error=$(test_errs[end]) convergence=$convergence")
        # models[i] = model
        error_string = inference_config.do_cross_validation ?
            "$alpha\t$beta\t$mu\t$v\t$(inference_config.fold_num)\t$(errs[end])\t$(test_errs[end])\t$(lls[end])\n" :
            "$alpha\t$beta\t$mu\t$v\t$(errs[end])\t$(test_errs[end])\t$(lls[end])\n"
        all_errors[i] = error_string

        (p, r) = predict_outcomes(pm, params, collect(values(pm.data.cell_lines)))
        log_message("done predicting")
        # predictions[i] = p
        # ranks[i] = r

        generic_filename =  "alpha_$(alpha)_beta_$(beta)_mean_$(mu)_var_$(v).txt"
        # writing result files
        write_results(pm, params, result_dir, generic_filename, p, r)

        resulting_ranking = joinpath(result_dir, "ranking", generic_filename)

        if inference_config.compute_wpc_index
            wpc_script_cmd = `perl $script_name $sd_file $test_response_file $zscore_file $resulting_ranking $wpc_target_file $wpc_summary_file`
            run(wpc_script_cmd)
        end

        p = nothing
        r = nothing
        model = nothing
        lls = nothing
        errs = nothing
        test_errs = nothing

        Base.GC.gc(true)
    end

    # write the errors file
    f = open(joinpath(inference_config.target_dir, "errors.txt"), "a")
    for s in all_errors
        @printf(f, "%s", s)
    end
    close(f)

end


text:
export  Scalar1D
mutable struct Scalar1D
    κ::Float64
end

function Scalar1D(prop::Dict{String, Any})
    κ = prop["kappa"]
    Scalar1D(κ)
end

function getStress(mat::Scalar1D, ∇u::Array{Float64}, ∇u_t::Array{Float64}, Δt::Float64 = 0.0)
    sigma = mat.κ * ∇u
    H = mat.κ * diagm(0=>ones(2))
    return sigma, H
end


text:
module Surrogates
using LinearAlgebra
using GLM
using Distributions
using Sobol
using LatinHypercubeSampling
using Requires

abstract type AbstractSurrogate <: Function end
include("Radials.jl")
include("Kriging.jl")
include("Sampling.jl")
include("Optimization.jl")
include("Lobachesky.jl")
include("LinearSurrogate.jl")
include("InverseDistanceSurrogate.jl")
include("SecondOrderPolynomialSurrogate.jl")

remove_tracker(x) = x

function __init__()
    @require XGBoost="009559a3-9522-5dbb-924b-0b6ed2b22bb9" begin
        using XGBoost
        include("RandomForestSurrogate.jl")
    end

    @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" begin
        using Flux
        using Flux: @epochs
        include("NeuralSurrogate.jl")
    end

    @require Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
        remove_tracker(x::TrackedReal) = Tracker.data(x)
        remove_tracker(x::TrackedArray) = Tracker.data(x)
    end

    @require LIBSVM="b1bec4e5-fd48-53fe-b0cb-9723c09d164b" begin
        using LIBSVM
        include("SVMSurrogate.jl")
    end
end

export AbstractSurrogate, SamplingAlgorithm
export Kriging, RadialBasis, add_point!, current_estimate, std_error_at_point
export sample, GridSample, UniformSample, SobolSample, LatinHypercubeSample, LowDiscrepancySample
export RandomSample
export SRBF,LCBS,EI,DYCORS,SOP,surrogate_optimize
export LobacheskySurrogate, lobachesky_integral, lobachesky_integrate_dimension
export LinearSurrogate
export RandomForestSurrogate
export SVMSurrogate
export NeuralSurrogate
export InverseDistanceSurrogate
export SecondOrderPolynomialSurrogate
end


text:
module SimulatedAnnealing

using Parameters
using Statistics

export
    AnnealingOptimization, AnnealingState,
    OVGCriterion, SSVCriterion,
    ConstantDecrementRule, AVLDecrementRule,
    simulated_annealing

"""
    AnnealingOptimization

Object representing an optimization done using a Simulated Annealing
algorithm.

Implement single stage homogeneous Simulated Annealing as (very well) described
in
Varanelli, 1996, *On the acceleration of simulated annealing*, Chapter 2.

Notation and naming convention is consistent with that work.

A configuration of the system being optimized is represented by the parametric
type `C`. The configurations to optimize can be of any type.

The object is iterable to allow access to the internal state during the
optimization process.

Parameters
==========
energy: callable(::C)::T
    Callable returning the energy of a configuration.
propose_candiate: callable(::C)::Tuple{::C, ::T}
    Callable returning a new configuration together based on an existing one
    togetehr with the energy difference between the two.
stop_criterion: callable(::AnnealingState)::Bool
    Callable returning whether the annealing must stop.
    Default is Otten-van Ginneken adaptive stop criterion.
decrement_rule: callable(::AneealingState)::T
    Callable returning the new temperature to use at the end of a Markov step.
    Default is Aarts and van Laarhoven decrement rule.
initial_temperature: Real
initial_configuration: C
neighborhood_size: integer
    Number of states that can be accessed from a given one.
"""
struct AnnealingOptimization{T, C, EFUNC, PFUNC, SFUNC, TFUNC}
    energy::EFUNC
    propose_candidate::PFUNC
    stop_criterion::SFUNC
    decrement_rule::TFUNC
    initial_temperature::T
    initial_configuration::C
    neighborhood_size::Int
end

function AnnealingOptimization(
        energy,
        propose_candidate,
        initial_configuration,
        initial_temperature::Real,
        neighborhood_size::Integer)

    return AnnealingOptimization(
        energy,
        propose_candidate,
        SSVCriterion(),
        AVLDecrementRule(),
        initial_temperature,
        initial_configuration,
        neighborhood_size)
end

"""
    AnnealingOptimization(energy, propose_candidate, samples, neighborhood_size)

Given samples of the configuration space, determine the initial temperature.
The first sample is used as initial configuration.

Initial temperature is the standard deviation of the energy in the sample,
according to White criterion (see eq. 2.36 in Varanelli).
"""
function AnnealingOptimization(energy, propose_candidate, samples::Vector,
                               neighborhood_size::Integer)

    energies = energy.(samples)
    return AnnealingOptimization(
        energy,
        propose_candidate,
        samples[1],
        std(energies),
        neighborhood_size)
end

"""
    AnnealingState

Struct representing the sate of the Simulated Annealing algorithm at the end
of a Markov chain step.

Fields
======
temperature
current_configuration
current_energy
bsf_configuration
    Best configuration found so far.
bsf_energy
    Energy of the best configuratin found so far.
energies
    Array of all energies encountered during the last Markov chain step.
"""
struct AnnealingState{T, C}
    temperature::T
    current_configuration::C
    current_energy::T
    bsf_configuration::C
    bsf_energy::T
    energies::Vector{T}
end

function initial_state(search::AnnealingOptimization{T}) where T
    config = search.initial_configuration
    E = search.energy(config)
    return AnnealingState(search.initial_temperature,
                          config, E,
                          config, E,
                          zeros(T, 0))
end

include("decrement_rule.jl")
include("stop_criterion.jl")

Base.eltype(::Type{A}) where {T, C, A <: AnnealingOptimization{T, C}} = AnnealingState{T, C}
Base.IteratorSize(::Type{A}) where {A <: AnnealingOptimization} = Base.SizeUnknown()

function Base.iterate(
        search::AnnealingOptimization{T},
        state=initial_state(search)) where T
    length(state.energies) > 0 && search.stop_criterion(state) && return nothing

    # Initialize all internal loop variables
    configuration = state.current_configuration
    E = state.current_energy
    bsf = state.bsf_configuration
    bsf_energy = state.bsf_energy
    energies = zeros(T, search.neighborhood_size)

    # Markov chain at constant temperature
    for k in 1:search.neighborhood_size
        candidate, dE = search.propose_candidate(configuration)

        if dE < 0 || rand() < exp(-dE/state.temperature)
            configuration = candidate
            E += dE

            if E < state.bsf_energy
                bsf = configuration
                bsf_energy = E
            end
        end

        energies[k] = E
    end

    new_state = AnnealingState(search.decrement_rule(state),
                               configuration, E,
                               bsf, bsf_energy,
                               energies)

    return new_state, new_state
end

function simulated_annealing(search::AnnealingOptimization)
    local state = nothing  # Avoid being shadowed inside the loop

    for s in search  # Go to the end of the search
        state = s
    end

    return state.bsf_configuration, state.bsf_energy
end

end


text:
# patches are characterized the by following four parameters
type Patch
    β::Float64 # MC's probability of predation
    λ::Float64 # MC's probability of finding food
    α::Int64 # MC's cost to visit
    Y::Int64 # MC's value of food
end

# an Ecosystem is essentially an array of patches
typealias Ecosystem Array{Patch,1};

# the patches from page 54
myplace = Ecosystem([
    Patch(0., 0., 1, 0),
    Patch(0.004, 0.4, 1, 3),
    Patch(0.02, 0.6, 1, 5)
    ]);

x_c = 3::Int64; # critical energy threshold
C = 10::Int64; # max capacity
T = 20::Int64; # time horizon

function xprime(x::Int64, α::Int64, y::Int64, x_critical::Int64, C::Int64)
    return clamp(x - α + y, x_critical, C)
end;

function xprimeprime(x::Int64, α::Int64, x_critical::Int64, C::Int64)
    return clamp(x - α, x_critical, C)
end;

function findOptimalSelections(
    aplace::Ecosystem,
    C::Int64,
    T::Int64,
    x_c::Int64)
    Fitness = zeros(Float64, C+1, T)
    Selection = zeros(Int64, C+1, T)
    for t = T:-1:1
        for x = 0:C
            if t == T
                Fitness[x+1, t] = x > x_c ? 1.0 : 0.0
            elseif x <= x_c
                Fitness[x+1, t] = 0.0
            else
                patchoptions = [(1-P.β)*
                    (P.λ*Fitness[xprime(x, P.α, P.Y, x_c, C)+1,t+1] +
                    (1.0-P.λ)*Fitness[xprimeprime(x, P.α, x_c, C)+1,t+1])
                    for P in aplace]
                (Fitness[x+1, t], Selection[x+1, t]) = findmax(patchoptions)
            end
        end
    end
    return Selection, Fitness
end;

(Selection, Fitness) = findOptimalSelections(myplace, C, T, x_c);

## compare to Mangel and Clark page 55
Fitness[5:11,19:-1:1]


text:
addChildrenVariable!(x::IntVar, y::IntVarView) = push!(x.children, y)
addChildrenVariable!(x::IntVarView, y::IntVarView) = addChildrenVariable!(x.x, y)

struct IntDomainViewMul <: IntDomainView
    orig            ::AbstractIntDomain
    a               ::Int
end

struct IntVarViewMul <: IntVarView
    x               ::AbstractIntVar
    a               ::Int
    domain          ::IntDomainViewMul
    id              ::String

    """
        IntVarViewMul(x::AbstractIntVar, a::Int, id::String)

    Create a *fake* variable `y`, such that `y == a*x`. This variable behaves like an usual one.
    """
    function IntVarViewMul(x::AbstractIntVar, a::Int, id::String)
        @assert a > 0
        dom = IntDomainViewMul(x.domain, a)
        var = new(x, a, dom, id)
        addChildrenVariable!(x, var)
        return var
    end
end

struct IntDomainViewOpposite <: IntDomainView
    orig            ::AbstractIntDomain
end

struct IntVarViewOpposite <: IntVarView
    x               ::AbstractIntVar
    domain          ::IntDomainViewOpposite
    id              ::String

    """
    IntVarViewOpposite(x::AbstractIntVar, id::String)

    Create a *fake* variable `y`, such that `y = -x`. This variable behaves like an usual one.
    """
    function IntVarViewOpposite(x::AbstractIntVar, id::String)
        dom = IntDomainViewOpposite(x.domain)
        var = new(x, dom, id)
        addChildrenVariable!(x, var)
        return var
    end
end

struct IntDomainViewOffset <: IntDomainView
    orig            ::AbstractIntDomain
    c               ::Int
end

struct IntVarViewOffset <: IntVarView
    x               ::AbstractIntVar
    c               ::Int
    domain          ::IntDomainViewOffset
    id              ::String

    """
    IntVarViewOffset(x::AbstractIntVar, id::String)

    Create a *fake* variable `y`, such that `y = x + c`. This variable behaves like an usual one.
    """
    function IntVarViewOffset(x::AbstractIntVar, c::Int, id::String)
        dom = IntDomainViewOffset(x.domain, c)
        var = new(x, c, dom, id)
        addChildrenVariable!(x, var)
        return var
    end
end

assignedValue(x::IntVarViewMul) = x.a * assignedValue(x.x)
assignedValue(x::IntVarViewOpposite) = -1 * assignedValue(x.x)
assignedValue(x::IntVarViewOffset) = assignedValue(x.x) + x.c

"""
    isempty(dom::IntDomainView)

Return `true` iff `dom` is an empty set.
"""
Base.isempty(dom::IntDomainView) = isempty(dom.orig)

reset_domain!(dom::IntDomainView) = reset_domain!(dom.orig)

"""
    rootVariable(x::IntVarView)

Return the "true" variable behind `x`.
"""
rootVariable(x::IntVarView) = rootVariable(x.x)

"""
    length(dom::IntDomainView)

Return the size of `dom`.
"""
Base.length(dom::IntDomainView) = length(dom.orig)

"""
    Base.in(value::Int, dom::IntDomainView)

Check if an integer is in the domain.
"""
function Base.in(value::Int, dom::IntDomainViewMul)
    if value % dom.a != 0
        return false
    end
    return (value ÷ dom.a) in dom.orig
end
Base.in(value::Int, dom::IntDomainViewOpposite) = -value in dom.orig
Base.in(value::Int, dom::IntDomainViewOffset) = value - dom.c in dom.orig

"""
    remove!(dom::IntDomainView, value::Int)

Remove `value` from `dom`.
"""
function remove!(dom::IntDomainViewMul, value::Int)
    if value % dom.a != 0
        return Int[]
    end

    return remove!(dom.orig, value ÷ dom.a) * dom.a
end
remove!(dom::IntDomainViewOpposite, value::Int) = -1 * remove!(dom.orig, -value)
remove!(dom::IntDomainViewOffset, value::Int) = map(x -> x + dom.c, remove!(dom.orig, value - dom.c))

"""
    removeAll!(dom::IntDomainView)

Remove every value from `dom`. Return the removed values.
"""
removeAll!(dom::IntDomainViewMul) = dom.a * removeAll!(dom.orig)
removeAll!(dom::IntDomainViewOpposite) = -1 * removeAll!(dom.orig)
removeAll!(dom::IntDomainViewOffset) = map(x -> x + dom.c, removeAll!(dom.orig))

"""
    minimum(dom::IntDomainView)

Return the minimum value of `dom`.
"""
minimum(dom::IntDomainViewMul) = dom.a * minimum(dom.orig)
minimum(dom::IntDomainViewOpposite) = -1 * maximum(dom.orig)
minimum(dom::IntDomainViewOffset) = dom.c + minimum(dom.orig)

"""
    maximum(dom::IntDomainView)

Return the maximum value of `dom`.
"""
maximum(dom::IntDomainViewMul) = dom.a * maximum(dom.orig)
maximum(dom::IntDomainViewOpposite) = -1 * minimum(dom.orig)
maximum(dom::IntDomainViewOffset) = dom.c + maximum(dom.orig)

"""
    removeAbove!(dom::IntDomainView, value::Int)

Remove every integer of `dom` that is *strictly* above `value`.
"""
removeAbove!(dom::IntDomainViewMul, value::Int) = dom.a * removeAbove!(dom.orig, convert(Int, floor(value / dom.a)))
removeAbove!(dom::IntDomainViewOpposite, value::Int) = -1 * removeBelow!(dom.orig, -value)
removeAbove!(dom::IntDomainViewOffset, value::Int) = map(x -> x + dom.c, removeAbove!(dom.orig, value - dom.c))

"""
    removeBelow!(dom::IntDomainView, value::Int)

Remove every integer of `dom` that is *strictly* below `value`. Return the pruned values.
"""
removeBelow!(dom::IntDomainViewMul, value::Int) = dom.a * removeBelow!(dom.orig, convert(Int, ceil(value / dom.a)))
removeBelow!(dom::IntDomainViewOpposite, value::Int) = -1 * removeAbove!(dom.orig, -value)
removeBelow!(dom::IntDomainViewOffset, value::Int) = map(x -> x + dom.c, removeBelow!(dom.orig, value - dom.c))

"""
    assign!(dom::IntDomainView, value::Int)

Remove everything from the domain but `value`. Return the removed values. Return the pruned values.
"""
function assign!(dom::IntDomainViewMul, value::Int)
    @assert value % dom.a == 0
    return dom.a * assign!(dom.orig, value ÷ dom.a)
end
assign!(dom::IntDomainViewOpposite, value::Int) = -1 * assign!(dom.orig, -value)
assign!(dom::IntDomainViewOffset, value::Int) = map(x -> x + dom.c, assign!(dom.orig, value - dom.c))

"""
    Base.iterate(dom::IntDomainView, state=1)

Iterate over the domain in an efficient way. The order may not be consistent.
WARNING: Do **NOT** update the domain you are iterating on.
"""
function Base.iterate(dom::IntDomainViewMul, state=1)
    returned = iterate(dom.orig, state)
    if isnothing(returned)
        return nothing
    end

    value, newState = returned
    return value * dom.a, newState
end

function Base.iterate(dom::IntDomainViewOpposite, state=1)
    returned = iterate(dom.orig, state)
    if isnothing(returned)
        return nothing
    end

    value, newState = returned
    return -value, newState
end

function Base.iterate(dom::IntDomainViewOffset, state=1)
    returned = iterate(dom.orig, state)
    if isnothing(returned)
        return nothing
    end

    value, newState = returned
    return value + dom.c, newState
end

"""
    updateMaxFromRemovedVal!(dom::IntDomainView, v::Int)

Knowing that `v` just got removed from `dom`, update `dom`'s maximum value.
"""
function updateMaxFromRemovedVal!(dom::IntDomainViewMul, v::Int)
    if maximum(dom) == v
        updateMaxFromRemovedVal!(dom.orig, v ÷ dom.a)
    end
end

function updateMaxFromRemovedVal!(dom::IntDomainViewOpposite, v::Int)
    if maximum(dom) == v
        updateMinFromRemovedVal!(dom.orig, -v)
    end
end

function updateMaxFromRemovedVal!(dom::IntDomainViewOffset, v::Int)
    if maximum(dom) == v
        updateMaxFromRemovedVal!(dom.orig, v - dom.c)
    end
end

"""
    updateMinFromRemovedVal!(dom::IntDomainView, v::Int)

Knowing that `v` just got removed from `dom`, update `dom`'s minimum value.
"""
function updateMinFromRemovedVal!(dom::IntDomainViewMul, v::Int)
    if minimum(dom) == v
        updateMinFromRemovedVal!(dom.orig, v ÷ dom.a)
    end
end

function updateMinFromRemovedVal!(dom::IntDomainViewOpposite, v::Int)
    if minimum(dom) == v
        updateMaxFromRemovedVal!(dom.orig, -v)
    end
end

function updateMinFromRemovedVal!(dom::IntDomainViewOffset, v::Int)
    if minimum(dom) == v
        updateMinFromRemovedVal!(dom.orig, v - dom.c)
    end
end

parentValue(y::IntVarViewMul, v::Int) = v ÷ y.a
parentValue(y::IntVarViewOffset, v::Int) = v - y.c
parentValue(::IntVarViewOpposite, v::Int) = - v
childrenValue(::IntVar, v::Int) = v
childrenValue(y::IntVarViewMul, v::Int) = childrenValue(y.x, v) * y.a
childrenValue(y::IntVarViewOffset, v::Int) = childrenValue(y.x, v) + y.c
childrenValue(y::IntVarViewOpposite, v::Int) = - childrenValue(y.x, v)


text:
using Query
using DataFrames

df = DataFrame(name=["John", "Sally", "Kirk"], age=[23., 42., 59.], children=[3,5,2])

x = @from i in df begin
    @orderby i.age
    @select {Name=lowercase(i.name)}
    @collect DataFrame
end

println(x)

x = @from i in df begin
    @orderby descending(i.age)
    @select {Name=lowercase(i.name)}
    @collect DataFrame
end

println(x)

x = @from i in df begin
    @orderby ascending(i.age)
    @select {Name=lowercase(i.name)}
    @collect DataFrame
end

println(x)


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ c2d0d93c-8d12-4b58-852b-10795ff1f8db
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ c74a1698-83dd-11eb-1bd2-c145ea201500
begin
	using DataFrames
	using LogisticRegression

	using Plots
	using StatsPlots
	gr()
end;

# ╔═╡ 43f669b2-8596-11eb-2e12-efdef4ac0e29
md"
# Demo of Logistic Regression

JingYu Nng
"

# ╔═╡ 135d298a-83de-11eb-1f50-83e4261f2bdf
md"
### Generate mock data
"

# ╔═╡ f374b00a-83dd-11eb-0f93-612e3ed66e75
begin
	n = 100
    β = [1., 1.]
    df = DataFrame(
        X₁=vcat(randn(n).-5, randn(n).+10),
        X₂=vcat(randn(n).-8, randn(n).+16),
        Y=vcat(ones(n), zeros(n))
    )
end

# ╔═╡ 2b9fc1a6-83de-11eb-3863-19df8f1b6252
md"
### Initial model
"

# ╔═╡ 25ae7788-83de-11eb-27be-2d51a3c90c8c
begin
	lrm = LogisticRegressionModel(df, :Y, [:X₁, :X₂])
end;

# ╔═╡ 7a44b194-83df-11eb-2de4-7f16b19d1239
begin
	const C_GRAD = cgrad([
		RGBA(53/255, 157/255, 219/255, 1),
		RGBA(240/255, 240/255, 240/255, 1),
		RGBA(219/255, 64/255, 68/255, 1)
	])

	function plot_model(model::LogisticRegressionModel, df::DataFrame; margin=3)
		lim = maximum(abs.(lrm.xs))
		x = -lim-margin:0.1:lim+margin
		heatmap(
			x, x, (x1, x2)->predict(lrm, [x1, x2]),
			title="Logistic Regression Model",
			xlabel="X₁",
			ylabel="X₂",
			legend=:topleft,
			clim=(0, 1),
			c=C_GRAD,
		)
		@df df[df.Y .== 0, :] scatter!(:X₁ , :X₂, color=C_GRAD[1], label="False")
		@df df[df.Y .== 1, :] scatter!(:X₁ , :X₂, color=C_GRAD[3], label="True")
	end
end;

# ╔═╡ da6bb3a2-83de-11eb-3a43-55dbeaabac48
md"
### Model before training
"

# ╔═╡ e78abf10-83de-11eb-29ca-c72ce80b5877
plot_model(lrm, df)

# ╔═╡ 31343290-83df-11eb-0c9f-114c43fd1c8d
md"
### Training
"

# ╔═╡ 37c5a346-83df-11eb-06f8-9b549328b45f
fit!(lrm, η=1e-2, atol=-1e-4);

# ╔═╡ 4242583c-83df-11eb-08b7-752569c5bbe8
md"
### Model after training
"

# ╔═╡ 8a450b0c-83df-11eb-34b9-4351fb9d772c
plot_model(lrm, df)

# ╔═╡ Cell order:
# ╟─43f669b2-8596-11eb-2e12-efdef4ac0e29
# ╟─c2d0d93c-8d12-4b58-852b-10795ff1f8db
# ╠═c74a1698-83dd-11eb-1bd2-c145ea201500
# ╟─135d298a-83de-11eb-1f50-83e4261f2bdf
# ╠═f374b00a-83dd-11eb-0f93-612e3ed66e75
# ╟─2b9fc1a6-83de-11eb-3863-19df8f1b6252
# ╠═25ae7788-83de-11eb-27be-2d51a3c90c8c
# ╟─7a44b194-83df-11eb-2de4-7f16b19d1239
# ╟─da6bb3a2-83de-11eb-3a43-55dbeaabac48
# ╠═e78abf10-83de-11eb-29ca-c72ce80b5877
# ╟─31343290-83df-11eb-0c9f-114c43fd1c8d
# ╠═37c5a346-83df-11eb-06f8-9b549328b45f
# ╟─4242583c-83df-11eb-08b7-752569c5bbe8
# ╠═8a450b0c-83df-11eb-34b9-4351fb9d772c


text:
module Player

abstract type Player end

function choose (current :: Player) :: Int
end

end

max_stars_repo_name: TheGluck/Connect-Four-AI

max_issues_repo_name: TheGluck/Connect-Four-AI

max_forks_repo_name: TheGluck/Connect-Four-AI


text:
# included by src.jl
# core types used in html constructing

# ------------------- type tree
abstract type AbstractHtmlTag <: Any end

# ------------------- public constants
const CONS = (
    redasterisk = "*", # red asterisk for required terms
)

# ------------------- synonym def
StrOrTag = Union{String,AbstractHtmlTag}
StrOrNo = Union{String,Nothing}
HtmlTagVector = Vector{T} where T <: AbstractHtmlTag

# ------------------- 成对tag: <tag>content</tag>
mutable struct PairedHtmlTag <: AbstractHtmlTag
    tag::String  # tag name, e.g. table, head
    id::StrOrNo  # tag id
    class::StrOrNo  # tag class applied
    style::StrOrNo  # CSS style
    singleAttributes::Vector{String}  # single attribute, e.g. required
    pairedAttributes::Dict{String,String}  # paired attributes, e.g. min=2000, width=\"100\"
    content::HtmlTagVector  # content between <tag> and </tag>
    # -------------
    function PairedHtmlTag(tag::String, id::StrOrNo, class::StrOrNo, style::StrOrNo;
        singleAttributes::Vector{String} = String[],
        pairedAttributes::Dict{String,String} = Dict{String,String}(),
        content::HtmlTagVector = AbstractHtmlTag[] )
        tag == "" ?  (throw(ErrorException("empty tag name found"))) : nothing
        new(tag,id,class,style,singleAttributes,pairedAttributes,content)
    end
end
# ------------------- 单个标签tag: <tag />
mutable struct SingleHtmlTag <: AbstractHtmlTag
    tag::String  # tag name, e.g. table, head
    id::StrOrNo  # tag id
    class::StrOrNo  # tag class applied
    style::StrOrNo  # CSS style
    singleAttributes::Vector{String}  # single attribute, e.g. required
    pairedAttributes::Dict{String,String}  # paired attributes, e.g. min=2000, width=\"100\"
    # -------------
    function SingleHtmlTag(tag::String, id::StrOrNo, class::StrOrNo, style::StrOrNo;
        singleAttributes::Vector{String} = String[],
        pairedAttributes::Dict{String,String} = Dict{String,String}() )
        tag == "" ?  (throw(ErrorException("empty tag name found"))) : nothing
        new(tag,id,class,style,singleAttributes,pairedAttributes)
    end
end
# -------------------- 空白tag：blank tag, (a placeolder for pure string contents)
struct BlankHtmlTag <: AbstractHtmlTag
    content::String
end
# -------------------- HTML page (.html file)
"""
    HtmlPage

type to define an HTML page file.
"""
mutable struct HtmlPage <: Any
    filename::String # file name when saved
    pagetitle::String # page title displayed on browser's status bar
    lang::String # language, "en-us" by default
    charset::String # <meta charset = "utf-8" /> by default
    head::Vector # a list of the elements in the  part
    body::Vector # a list of the elements in the  part, ordered
    # ---------------------------
    function HtmlPage( filename::String = "index.html",
        pagetitle::String = "main", lang::String = "en-us",
        charset::String = "utf-8" )
        new(filename,pagetitle,lang,charset,[],[])
    end
end

#

max_stars_repo_name: Clpr/HtmlConstructor

max_issues_repo_name: Clpr/HtmlConstructor

max_forks_repo_name: Clpr/HtmlConstructor


text:
#
# newuoa.jl --
#
# Julia interface to Mike Powell's NEWUOA method.
#
# -----------------------------------------------------------------------------
#
# This file is part of OptimPack.jl which is licensed under the MIT
# "Expat" License:
#
# Copyright (C) 2015-2019, Éric Thiébaut <https://github.com/emmt/OptimPack.jl>.
#

module Newuoa

export
    newuoa,
    newuoa!

using Compat
using Compat.Printf

import
    ..AbstractContext,
    ..AbstractStatus,
    ..getncalls,
    ..getradius,
    ..getreason,
    ..getstatus,
    ..grow!,
    ..iterate,
    ..restart

# The dynamic library implementing the method.
import ..libnewuoa

# Status returned by most functions of the library.
struct Status <: AbstractStatus
    _code::Cint
end

# Possible status values returned by NEWUOA.
const INITIAL_ITERATE      = Status( 2)
const ITERATE              = Status( 1)
const SUCCESS              = Status( 0)
const BAD_NVARS            = Status(-1)
const BAD_NPT              = Status(-2)
const BAD_RHO_RANGE        = Status(-3)
const BAD_SCALING          = Status(-4)
const ROUNDING_ERRORS      = Status(-5)
const TOO_MANY_EVALUATIONS = Status(-6)
const STEP_FAILED          = Status(-7)
const BAD_ADDRESS          = Status(-8)
const CORRUPTED            = Status(-9)

# Get a textual explanation of the status returned by NEWUOA.
function getreason(status::Status)
    ptr = ccall((:newuoa_reason, libnewuoa), Ptr{UInt8}, (Cint,), status._code)
    if ptr == C_NULL
        error("unknown NEWUOA status: ", status._code)
    end
    unsafe_string(ptr)
end

"""
# Minimizing a function of many variables

Mike Powell's **NEWUOA** algorithm attempts to find the variables `x` which
solve the problem:

    min f(x)

where `x` is a vector of variables that has `n ≥ 2` components and `f(x)` is
an objective function.  The algorithm employs quadratic approximations to the
objective which interpolates the objective function at `m` points, the value
`m = 2n + 1` being recommended.  The parameter `rho` controls the size of the
trust region and it is reduced automatically from `rhobeg` to `rhoend` (such
that `0 < rhoend ≤ rhobeg`).

The in-place version of the algorithm is called as:

    Newuoa.minimize!(f, x, rhobeg, rhoend) -> (status, x, fx)

where `f` is the objective function, `x` is a vector with the initial and final
variables, `rhobeg` and `rhoend` are the initial and final sizes of the trust
region.  The result is a tuple of 3 values: `status` indicates whether the
algorithm was successful, `x` is the final value of the variables and `fx =
f(x)` is the objective function at `x`.  Normally, `status` should be
`Newuoa.SUCCESS`; otherwise, `getreason(status)` yields a textual explanation
about the failure.

The method:

    Newuoa.minimize(f, x0, rhobeg, rhoend) -> (status, x, fx)

is identical to `Newuoa.minimize!` but does not modify the vector `x0` of
initial variables.

## Precision and scaling of variables

Parameter `rhobeg` should be set to the typical size (in terms of Euclidean
norm of the change of variables) of the region to explorate and `rhoend`
should be set to the typical precision.  The proper scaling of the variables is
important for the success of the algorithm and the optional `scale` keyword
should be specified if the typical precision is not the same for all variables.
If specified, `scale` is an array of strictly nonnegative values and of same
size as the variables `x`, such that `scale[i]*rho` (with `rho` the trust
region radius) is the size of the trust region for the `i`-th variable.  If
keyword `scale` is not specified, a unit scaling for all the variables is
assumed.

## Keywords

The following keywords are available:

* `scale` specifies the typical magnitudes of the variables.  If specified, it
  must have as many elements as `x`, all strictly positive.

* `check` (`true` by default) specifies whether to throw an exception if the
  algorithm is not fully successful.

* `verbose` (`0` by default) set the amount of printing.

* `maxeval` specifies the maximum number of calls to the objective function.
  The default setting is `maxeval = 30n` with `n = length(x)` the number of
  variables.

* `npt` specifies the number of points to use for the quadratic approximation
  of the objective function.  The default setting is the recommended value:
  `npt = 2n + 1` with `n = length(x)` the number of variables.

* `work` specifies a workspace to (re)use.  It must be a vector of double
  precision floating-point values.  If it is too small, its size is
  automatically adjusted (by calling [`resize!`](@ref)).  This keyword is
  useful to avoid any new allocation (and garbage colection) when several
  similar optimizations are to be performed.

## References

The algorithm is described in:

* M.J.D. Powell, "The NEWUOA software for unconstrained minimization without
  derivatives," in Large-Scale Nonlinear Optimization, editors G. Di Pillo and
  M. Roma, Springer, pp. 255-297 (2006).

"""
minimize(args...; kwds...) = optimize(args...; maximize=false, kwds...)
minimize!(args...; kwds...) = optimize!(args...; maximize=false, kwds...)
@doc @doc(minimize) minimize!

"""

    Newuoa.maximize(f, x0, rhobeg, rhoend) -> (status, x, fx)
    Newuoa.maximize!(f, x, rhobeg, rhoend) -> (status, x, fx)

are similar to `Newuoa.minimize` and `Newuoa.minimize!` respectively but
solve the unconstrained maximization problem:

    max f(x)

"""
maximize(args...; kwds...) = optimize(args...; maximize=true, kwds...)
maximize!(args...; kwds...) = optimize!(args...; maximize=true, kwds...)
@doc @doc(maximize) maximize!

# `_wrklen(...)` yields the number of elements in NEWUOA workspace.
_wrklen(n::Integer, npt::Integer) = _wrklen(Int(n), Int(npt))
_wrklen(n::Int, npt::Int) = (npt + 13)*(npt + n) + div(3*n*(n + 3),2)
_wrklen(x::AbstractVector{<:AbstractFloat}, npt::Integer) =
    _wrklen(length(x), npt)
function _wrklen(x::AbstractVector{<:AbstractFloat},
                 scl::AbstractVector{<:AbstractFloat},
                 npt::Integer)
    return _wrklen(x, npt) + length(scl)
end

# `_work(...)` yields a large enough workspace for NEWUOA.
_work(x::AbstractVector{<:AbstractFloat}, npt::Integer) =
    Vector{Cdouble}(undef, _wrklen(x, npt))
function _work(x::AbstractVector{<:AbstractFloat},
               scl::AbstractVector{<:AbstractFloat},
               npt::Integer)
    return Vector{Cdouble}(undef, _wrklen(x, scl, npt))
end

# Wrapper for the objective function in NEWUOA, the actual objective function
# is provided by the client data as a `jl_value_t*` pointer.
function _objfun(n::Cptrdiff_t, xptr::Ptr{Cdouble}, fptr::Ptr{Cvoid})::Cdouble
    x = unsafe_wrap(Array, xptr, n)
    f = unsafe_pointer_to_objref(fptr)
    return Cdouble(f(x))
end

# With precompilation, `__init__()` carries on initializations that must occur
# at runtime like `@cfunction` which returns a raw pointer.
const _objfun_c = Ref{Ptr{Cvoid}}()
function __init__()
    _objfun_c[] = @cfunction(_objfun, Cdouble,
                             (Cptrdiff_t, Ptr{Cdouble}, Ptr{Cvoid}))
end

"""
The methods:

    Newuoa.optimize(fc, x0, rhobeg, rhoend) -> (status, x, fx)
    Newuoa.optimize!(fc, x, rhobeg, rhoend) -> (status, x, fx)

are identical to `Newuoa.minimize` and `Newuoa.minimize!` respectively but have
an additional `maximize` keyword which is `false` by default and which
specifies whether to maximize the objective function; otherwise, the method
attempts to minimize the objective function.

"""
optimize(f::Function, x0::AbstractVector{<:Real}, args...; kwds...) =
    optimize!(f, copyto!(Array{Cdouble}(undef, length(x0)), x0),
              args...; kwds...)

function optimize!(f::Function, x::DenseVector{Cdouble},
                   rhobeg::Real, rhoend::Real;
                   scale::DenseVector{Cdouble} = Cdouble[],
                   maximize::Bool = false,
                   npt::Integer = 2*length(x) + 1,
                   check::Bool = true,
                   verbose::Integer = 0,
                   maxeval::Integer = 30*length(x),
                   work::Vector{Cdouble} = _work(x, scale, npt))
    n = length(x)
    nscl = length(scale)
    if nscl == 0
        sclptr = Ptr{Cdouble}(0)
    elseif nscl == n
        sclptr = pointer(scale)
    else
        error("bad number of scaling factors")
    end
    grow!(work, _wrklen(x, scale, npt))
    status = Status(ccall((:newuoa_optimize, libnewuoa), Cint,
                          (Cptrdiff_t, Cptrdiff_t, Cint, Ptr{Cvoid}, Any,
                           Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble,
                           Cptrdiff_t, Cptrdiff_t, Ptr{Cdouble}),
                          n, npt, maximize, _objfun_c[], f, x, sclptr,
                          rhobeg, rhoend, verbose, maxeval, work))
    if check && status != SUCCESS
        error(getreason(status))
    end
    return (status, x, work[1])
end

@doc @doc(optimize) optimize!

# Basic version similar to the FORTRAN version.
newuoa(f::Function, x0::DenseVector{Cdouble}, args...; kwds...) =
    newuoa!(f, copy(x0), args...; kwds...)

function newuoa!(f::Function, x::DenseVector{Cdouble},
                 rhobeg::Real, rhoend::Real;
                 npt::Integer = 2*length(x) + 1,
                 verbose::Integer = 0,
                 maxeval::Integer = 30*length(x),
                 check::Bool = true,
                 work::Vector{Cdouble} = _work(x, npt))
    n = length(x)
    grow!(work, _wrklen(x, npt))
    status = Status(ccall((:newuoa, libnewuoa), Cint,
                          (Cptrdiff_t, Cptrdiff_t, Ptr{Cvoid}, Any,
                           Ptr{Cdouble}, Cdouble, Cdouble, Cptrdiff_t,
                           Cptrdiff_t, Ptr{Cdouble}),
                          n, npt, _objfun_c[], f, x, rhobeg, rhoend,
                          verbose, maxeval, work))
    if check && status != SUCCESS
        error(getreason(status))
    end
    return (status, x, work[1])
end

"""

```julia
using OptimPack.Powell
ctx = Newuoa.create(n, rhobeg, rhoend; npt=..., verbose=..., maxeval=...)
```

creates a new reverse communication workspace for NEWUOA algorithm.  A typical
usage is:

```julia
x = Array{Cdouble}(undef, n)
x[...] = ... # initial solution
ctx = Newuoa.Context(n, rhobeg, rhoend; verbose=1, maxeval=500)
status = getstatus(ctx)
while status == Newuoa.ITERATE
    fx = ...       # compute function value at X
    status = iterate(ctx, fx, x)
end
if status != Newuoa.SUCCESS
    println("Something wrong occured in NEWUOA: ", getreason(status))
end
```

""" Context

# Context for reverse communication variant of NEWUOA.
# Must be mutable to be finalized.
mutable struct Context <: AbstractContext
    ptr::Ptr{Cvoid}
    n::Int
    npt::Int
    rhobeg::Cdouble
    rhoend::Cdouble
    verbose::Int
    maxeval::Int
    function Context(n::Integer, rhobeg::Real, rhoend::Real;
                     npt::Integer = 2*length(x) + 1,
                     verbose::Integer = 0,
                     maxeval::Integer = 30*length(x))
        ptr = ccall((:newuoa_create, libnewuoa), Ptr{Cvoid},
                    (Cptrdiff_t, Cptrdiff_t, Cdouble, Cdouble,
                     Cptrdiff_t, Cptrdiff_t),
                    n, npt, rhobeg, rhoend, verbose, maxeval)
        ptr != C_NULL || error(errno() == Base.Errno.ENOMEM
                               ? "insufficient memory"
                               : "invalid argument(s)")
        return finalizer(ctx -> ccall((:newuoa_delete, libnewuoa), Cvoid,
                                      (Ptr{Cvoid},), ctx.ptr),
                         new(ptr, n, npt, rhobeg, rhoend, verbose, maxeval))
    end
end

@deprecate create(args...; kwds...) Context(args...; kwds...)

function iterate(ctx::Context, f::Real, x::DenseVector{Cdouble})
    length(x) == ctx.n || error("bad number of variables")
    Status(ccall((:newuoa_iterate, libnewuoa), Cint,
                       (Ptr{Cvoid}, Cdouble, Ptr{Cdouble}),
                       ctx.ptr, f, x))
end

restart(ctx::Context) =
    Status(ccall((:newuoa_restart, libnewuoa), Cint, (Ptr{Cvoid},), ctx.ptr))

getstatus(ctx::Context) =
    Status(ccall((:newuoa_get_status, libnewuoa), Cint, (Ptr{Cvoid},), ctx.ptr))

getncalls(ctx::Context) =
    Int(ccall((:newuoa_get_nevals, libnewuoa), Cptrdiff_t, (Ptr{Cvoid},), ctx.ptr))

getradius(ctx::Context) =
    ccall((:newuoa_get_rho, libnewuoa), Cdouble, (Ptr{Cvoid},), ctx.ptr)

end # module Newuoa


text:
function getraster(T::Type{<:RasterDataSource}, layers=layers(T), args...; kw...)
    map(layers) do l
        _check_layer(T, l)
        getraster(T, l, args...; kw...)
    end
end

function _maybe_download(uri::URI, filepath)
    if !isfile(filepath)
        mkpath(dirname(filepath))
        println("Starting download for $uri")
        HTTP.download(string(uri), filepath)
    end
    filepath
end

function rasterpath()
    if haskey(ENV, "RASTERDATASOURCES_PATH") && isdir(ENV["RASTERDATASOURCES_PATH"])
        ENV["RASTERDATASOURCES_PATH"]
    else
        error("You must set `ENV[\"RASTERDATASOURCES_PATH\"]` to a path in your system")
    end
end

function delete_rasters()
    # May need an "are you sure"? - this could be a lot of GB of data to lose
    ispath(rasterpath()) && rm(rasterpath())
end

function delete_rasters(T::Type)
    ispath(rasterpath(T)) && rm(rasterpath(T))
end

_check_res(T, res) =
    res in resolutions(T) || throw(ArgumentError("Resolution $res not in $(resolutions(T))"))
_check_layer(T, layer) =
    layer in layers(T) || throw(ArgumentError("Layer $layer not in $(layers(T))"))

_date2string(t, date) = Dates.format(date, _dateformat(t))
_string2date(t, d::AbstractString) = Date(d, _dateformat(t))

_date_sequence(dates::AbstractArray, step) = dates
_date_sequence(dates::NTuple{2}, step) = first(dates):step:last(dates)
_date_sequence(date, step) = date:step:date


text:
using HTTP

"""
    HTTPStore

A basic HTTP store without any credentials. The underlying data is supposed to be
consolidated and only read operations are supported. This store is compatible to
datasets being served through the [xpublish](https://xpublish.readthedocs.io/en/latest/)
python package.
"""
struct HTTPStore <: AbstractStore
    url::String
end

function Base.getindex(s::HTTPStore, k::String)
r = HTTP.request("GET",string(s.url,"/",k),status_exception = false)
if r.status >= 300
    if r.status == 404
        nothing
    else
        error("Error connecting to $(s.url) :", String(r.body))
    end
else
    r.body
end
end

push!(storageregexlist,r"^https://"=>HTTPStore)
push!(storageregexlist,r"^http://"=>HTTPStore)
storefromstring(::Type{<:HTTPStore}, s,_) = ConsolidatedStore(HTTPStore(s),""),""

## This is a server implementation for Zarr datasets

function zarr_req_handler(s::AbstractStore, p)
  if s[p,".zmetadata"] === nothing
    consolidate_metadata(s)
  end
  request -> begin
    k = request.target
    k = lstrip(k,'/')
    contains("..",k) && return nothing
    r = s[p,k]
    try
      if r ===  nothing
        return HTTP.Response(404, "Error: Key $k not found")
      else
        return HTTP.Response(200, r)
      end
    catch e
      return HTTP.Response(404, "Error: $e")
    end
  end
end

HTTP.serve(s::AbstractStore, p, args...; kwargs...) = HTTP.serve(zarr_req_handler(s,p),args...;kwargs...)