text:
#!/usr/bin/env julia
using Libdl

dlopen("./libkashim", RTLD_GLOBAL)

struct Payload
    handle::Ptr{Cvoid}
    barrier::Ptr{Cvoid}
    function Payload(handle::Ptr{Cvoid})
        barrier = ccall(:ka_barrier_malloc, Ptr{Cvoid}, ())
        ccall(:ka_barrier_init, Cint, (Ptr{Cvoid}, Cuint), barrier, 2)
        new(handle, barrier)
    end
end

# Outside event the one we want to synchronize upon
cpu_event =  Base.Threads.Event()

# AsyncCondition is edge-triggered. We need to turn it into
# a level-triggered event. So we set up a barrier and have the AsyncCondition
# notify that, instead of waitingon the AsyncCondition directly.
barrier = Base.Threads.Event() # level-triggered
async = Base.AsyncCondition() do cond # edge-triggered
    notify(barrier)
    close(cond) # stop waiting for another notifycation
end
payload = Ref(Payload(async.handle))
@info "Setting up done"

# Task that is used to perform the waiting on cpu_event and the subsequent
# blocking barrier. We perhaps could use `@threadcall` to use LibUV threads
# to do the barrier_wait, but that barrier_wait should be cheap since the
# callback thread has already reached it.
T = Threads.@spawn begin
    GC.@preserve async payload begin
        # wait on MT-safe barrier
        @info "Waiting on notify"
        wait(barrier) # suspend ourselves until then
        @info "Received notify"

        # wait on actual work
        try
            @info "Waiting on cpu_event"
            wait(cpu_event)
        catch err
            bt = catch_backtrace()
            @error "Error thrown during ait on cpu_event" _ex=(err, bt)
        finally
            @info "Notifying callback"
            # now notify callback
            ka_barrier = payload[].barrier
            if ccall(:ka_barrier_wait, Cint, (Ptr{Cvoid},), ka_barrier) > 0
                ccall(:ka_barrier_destroy, Cvoid, (Ptr{Cvoid},), ka_barrier)
                ccall(:ka_barrier_free, Cvoid, (Ptr{Cvoid},), ka_barrier)
            end
            @info "done"
        end
    end
    return nothing
end

# This is the callback we want to synchronize with
# emulate a C-thread calling this by using `@threadcall` which
# executes a callback on the LibUV threadpool.
T2 = @async begin
    @info "threadcall callback"
    Base.@threadcall(:ka_callback, Cvoid, (Ptr{Cvoid},), payload)
    @info "callback done"
end

@info "notifying outside event"
notify(cpu_event)
@info "waiting upon T"
wait(T)
@info "waiting upon T2"
wait(T2)
@info "all done!"

max_stars_repo_name: vchuravy/KAShim

max_issues_repo_name: vchuravy/KAShim

max_forks_repo_name: vchuravy/KAShim


text:
################################################################
# List of headers to be automatically wrapped
# HEADER_FILES will be used by generator.jl
################################################################

const P4EST_HEADER_FILES = [
    "p4est.h",
    "p4est_bits.h",
    "p4est_connectivity.h",
    "p4est_extended.h",
    "p4est_geometry.h",
    "p4est_ghost.h",
    "p4est_iterate.h",
    "p4est_lnodes.h",
    "p4est_mesh.h",
    "p4est_search.h",
    "p4est_vtk.h"]

const P6EST_HEADER_FILES = [
    "p6est.h",
    "p6est_extended.h",
    "p6est_ghost.h"]

const P8EST_HEADER_FILES = [
    "p8est.h",
    "p8est_bits.h",
    "p8est_connectivity.h",
    "p8est_extended.h",
    "p8est_geometry.h",
    "p8est_ghost.h",
    "p8est_iterate.h",
    "p8est_lnodes.h",
    "p8est_mesh.h",
    "p8est_search.h",
    "p8est_vtk.h"]

const PXEST_HEADER_FILES = vcat(P4EST_HEADER_FILES,P6EST_HEADER_FILES,P8EST_HEADER_FILES)

const SC_HEADER_FILES = [
    "sc.h",
    "sc_containers.h",
    "sc_io.h"]

ext: jl


text:
abstract type Node end

struct Group <: Node
    root::String
    rel_path::String
    children::Vector{Node}
    config::Dict{Symbol,Any}
    function Group(root, rel_path, children, config)
        d = joinpath(root, rel_path)
        isdir(d) || throw(ArgumentError("Not a directory: $d"))
        check_config(config)
        new(root, rel_path, children, config)
    end
end

struct Document <: Node
    root::String
    rel_path::String
    kind::Symbol
    config::Dict{Symbol,Any}
    function Document(root, rel_path, children, config)
        f = joinpath(root, rel_path)
        isfile(f) || throw(ArgumentError("Not a file: $f"))
        check_config(config)
        new(root, rel_path, children, config)
    end
end

struct File <: Node
    root::String
    rel_path::String
end

const BUILDS = (:markdown, :script, :notebook)

const CONFIG_DEFAULTS = Dict(
    :active => true,
    :builds => BUILDS,
    :hide => false,
    :short_title => :use_title,
    :title => nothing,
    :weight => nothing,
    :deps => []
)

function group(root; skipliterate::Bool = false)
    config = parsefile_config(joinpath(root, "_config.jl"))
    inheritconfig!(config, CONFIG_DEFAULTS)

    grp = Group(root, ".", Node[], config)

    for child in readdir(root)
        if child == "_config.jl"
            continue
        elseif isdir(joinpath(root, child))
            push!(grp.children, group(child, grp; skipliterate = skipliterate))
        elseif endswith(child, ".jl") && !skipliterate
            push!(grp.children, document(child, grp))
        elseif endswith(child, ".md")
            push!(grp.children, document(child, grp))
        else
            push!(grp.children, file(child, grp))
        end
    end
    grp
end

function group(rel_path, parent::Group; skipliterate = false)
    config = parsefile_config(joinpath(parent.root, rel_path, "_config.jl"))
    inheritconfig!(config, parent.config)

    grp = Group(parent.root, rel_path, Node[], config)

    for child in readdir(joinpath(parent.root, rel_path))
        if child == "_config.jl"
            continue
        else
            child = joinpath(rel_path, child)
            if isdir(joinpath(parent.root, child))
                push!(grp.children, group(child, grp))
            elseif endswith(child, ".jl") && !skipliterate
                push!(grp.children, document(child, grp))
            elseif endswith(child, ".md")
                push!(grp.children, document(child, grp))
            else
                push!(grp.children, file(child, grp))
            end
        end
    end
    grp
end

function document(rel_path, parent::Group)
    kind, config, body = parsefile(joinpath(parent.root, rel_path))
    inheritconfig!(config, parent.config)
    Document(parent.root, rel_path, kind, config)
end

function file(rel_path, parent::Group)
    File(parent.root, rel_path)
end

function parsefile(path)
    try
        content = read(path, String)
        if isliterate(path, content)
            kind = :literate
            config, body = parse_literate(content)
        elseif isdocumenter(path, content)
            kind = :documenter
            config, body = parse_documenter(content)
        else
            error("Expected either a Documenter or Literate formatted file")
        end
        return kind, parse_config(config), body
    catch e
        error("""
             $(sprint(showerror, e))
             while processing $path
             """)
    end
end

function parse_documenter(content::String)
    md = Markdown.parse(content)
    body_blocks, content_blocks = Markdown.Code[], []
    for block in md.content
        if block isa Markdown.Code && block.language == "@cfg"
            push!(body_blocks, block)
        else
            push!(content_blocks, block)
        end
    end
    config = join(map(b -> b.code, body_blocks), '\n') * '\n'
    body = string(Markdown.MD(content_blocks))
    (config = config, body = body)
end

function parse_literate(content::String)
    lines = collect(eachline(IOBuffer(content)))
    config, body = String[], String[]
    for line in lines
        line = rstrip(line)
        m = match(r"^#cfg\h*(?<content>.+)", line)
        if m === nothing
            push!(body, line)
        else
            push!(config, rstrip(m[:content]))
        end
    end

    if isempty(config)
        error("No #cfg lines found")
    else
        config = join(config, '\n') * '\n'
        body = join(body, '\n') * '\n'
        return (config = config, body = body)
    end
end

function parsefile_config(path)
    try
        return parse_config(read(path, String))
    catch e
        error("""
             $(sprint(showerror, e))
             while parsing $path
             """)
    end
end

function parse_config(config_block::String)
    mod = execute_block(config_block)
    config = Dict{Symbol,Any}()
    for k in keys(CONFIG_DEFAULTS)
        isdefined(mod, k) && (config[k] = getfield(mod, k))
    end
    config
end

function inheritconfig!(dst, src)
    for k in keys(CONFIG_DEFAULTS)
        if k === :short_title
            dst[k] = get(dst, k, :use_title)
        else
            get(dst, k, nothing) === nothing && (dst[k] = src[k])
        end
    end
end

function check_config(config)
    c = config
    _check_type(k, T) = c[k] isa T || error("Config option `$k` must be of type `$T`")

    for k in keys(CONFIG_DEFAULTS)
        haskey(c, k) || error("Missing config option `$k`")
    end

    _check_type(:active, Bool)
    _check_type(:builds, Tuple)
    for b in c[:builds]
        b in BUILDS || error("Invalid build option `$b`. Valid options are: $BUILDS")
    end
    _check_type(:hide, Bool)
    _check_type(:title, AbstractString)
    if c[:short_title] === :use_title
        c[:short_title] = c[:title]
    elseif !(c[:short_title] isa AbstractString)
        error("Config option `$k` must be of type `AbstractString` or `:use_title`")
    end
    _check_type(:weight, Integer)
end

ismarkdown(path, content) = endswith(path, ".md")
isliterate(path, content) = endswith(path, ".jl")
isdocumenter(path, content) = ismarkdown(path, content)


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

name = "SoapyPlutoSDR"
version = v"0.2.1" # not yet tagged

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/pothosware/SoapyPlutoSDR.git", "ac9a9da5c14c73e752796618d56e259ca1ac6b11")
]

dependencies = [
    Dependency("libiio_jll"; compat="~0.23.0"),
    Dependency("soapysdr_jll"; compat="~0.8.0"),
    Dependency("libad9361_iio_jll"; compat="~0.2.0")
]

# Bash recipe for building across all platforms
script = raw"""
cd SoapyPlutoSDR

mkdir build && cd build
if [[ "${target}" != *-apple-* ]]; then
    cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
        -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
        -DCMAKE_BUILD_TYPE=Release \
        ..
    make -j${nproc}
    make install
fi
if [[ "${target}" == *-apple-* ]]; then
    cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
        -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
        -DCMAKE_BUILD_TYPE=Release \
        -DLibIIO_INCLUDE_DIR=${libdir}/iio.framework/Versions/0.23/Headers/ \
        -DLibIIO_LIBRARY=${libdir}/iio.framework/Versions/0.23/iio \
        -DLibAD9361_INCLUDE_DIR=${libdir}/ad9361.framework/Versions/0.2/Headers/ \
        -DLibAD9361_LIBRARY=${libdir}/ad9361.framework/Versions/0.2/ad9361 \
        ..
    make -j${nproc}
    make install
    mv ${libdir}/SoapySDR/modules0.8/libPlutoSDRSupport.so  ${libdir}/SoapySDR/modules0.8/libPlutoSDRSupport.dylib
fi
"""

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

# The products that we will ensure are always built
products = Product[
    LibraryProduct("libPlutoSDRSupport", :libPlutoSDRSupport, ["lib/SoapySDR/modules0.8/"])
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Example
using PackageWithBuildSpecificTestDeps
using Test

@test PackageWithBuildSpecificTestDeps.f(3) == 3

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
# ckl example
using Revise
using BridgeLandmarks
const BL = BridgeLandmarks
using Distributions
using LinearAlgebra
using Random
using FileIO

# using RCall
# using Random
# using DataFrames
# using DelimitedFiles
# using CSV
# using StaticArrays
#using JLD2

Random.seed!(9)

workdir = @__DIR__
cd(workdir)
outdir = workdir

#-------- read data ----------------------------------------------------------
dat = load("data_exp-ckl.jld2")
xobs0 = dat["xobs0"]
xobsT = dat["xobsT"]

p_ms = Pars_ms(δmom=0.1,cinit = 0.2, γinit = 2.0,ρinit = 0.95)
p_ahs = Pars_ahs(δmom=0.1,cinit = 0.02, γinit = 0.2, ρinit = 0.95, stdev=0.75, db=[2.5, 2.5])

ainit = mean(norm.(diff(xobs0)))
ups = [:innov, :mala_mom, :parameter]

@time landmarksmatching(xobs0,xobsT; ITER=250,pars=p_ms,outdir=outdir, ainit=ainit, updatescheme=ups)

plotlandmarksmatching(outdir)


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

""" Dynamics solver for JuliaFEM. """
module FEMDynamics

using Reexport
@reexport using FEMBase
using DifferentialEquations

type LinearImplicit <: AbstractAnalysis
    tspan :: Tuple{Float64,Float64}
    sol
    K :: SparseMatrixCSC{Float64}
    M :: SparseMatrixCSC{Float64}
    f :: SparseVector{Float64}
    u0 :: SparseVector{Float64}
    v0 :: SparseVector{Float64}
end

function LinearImplicit()
    tspan = (0.0, 1.0)
    return LinearImplicit(tspan, nothing, spzeros(0,0), spzeros(0,0),
                          spzeros(0), spzeros(0), spzeros(0))
end

"""
    get_global_matrices(analysis)

Assemble global matrices for problem.
"""
function get_global_matrices(analysis::Analysis{LinearImplicit})
    p = analysis.properties
    if !isempty(p.K) && !isempty(p.M)
        return p.K, p.M, p.f
    end

    time = first(analysis.properties.tspan)

    for problem in get_problems(analysis)
        assemble!(problem, time)
        is_field_problem(problem) && assemble_mass_matrix!(problem, time)
    end

    M = SparseMatrixCOO()
    K = SparseMatrixCOO()
    Kg = SparseMatrixCOO()
    f = SparseMatrixCOO()
    fg = SparseMatrixCOO()

    for problem in get_problems(analysis)
        is_field_problem(problem) || continue
        append!(M, problem.assembly.M)
        append!(K, problem.assembly.K)
        append!(Kg, problem.assembly.Kg)
        append!(f, problem.assembly.f)
        append!(fg, problem.assembly.fg)
    end

    dim = size(K, 1)

    M = sparse(M, dim, dim)
    K = sparse(K, dim, dim) + sparse(Kg, dim, dim)
    f = sparse(f, dim, 1) + sparse(fg, dim, 1)

    p.K = 1/2*(K + K')
    p.M = 1/2*(M + M')
    p.f = f
    return p.K, p.M, p.f
end

function FEMBase.run!(analysis::Analysis{LinearImplicit})

    # Preprocess: create global matrices, eliminate boundary conditions, remove
    # any nonzero rows and factorize

    K, M, f = get_global_matrices(analysis)
    ndofs = size(K, 1)
    full_vec = zeros(ndofs) # for temporary operations
    dim = get_unknown_field_dimension(first(get_problems(analysis)))
    nnodes = Int(ndofs/dim)

    for problem in get_problems(analysis)
        is_boundary_problem(problem) || continue
        eliminate_boundary_conditions!(problem, K, M, f)
    end
    dropzeros!(K)
    dropzeros!(M)

    if isempty(analysis.properties.u0)
        analysis.properties.u0 = spzeros(ndofs)
    end
    if isempty(analysis.properties.v0)
        analysis.properties.v0 = spzeros(ndofs)
    end
    u0 = analysis.properties.u0
    v0 = analysis.properties.v0

    nz = get_nonzero_rows(K)
    K = K[nz,nz]
    M = M[nz,nz]
    f = f[nz]
    u0 = u0[nz]
    v0 = v0[nz]
    M_fact = cholfact(M)
    ndofs = size(K, 1)

    # Solution

    function model(dx, x, p, t)
        info("Solving at time $t")
        ndofs = round(Int, length(dx)/2)
        u = x[1:ndofs]
        v = x[ndofs+1:end]
        dx[1:ndofs] = v
        dx[ndofs+1:end] = M_fact \ full(f - K*u)
    end

    x0 = [u0; v0]
    prob = ODEProblem(model, x0, analysis.properties.tspan)
    sol = analysis.properties.sol = solve(prob)

    # postprocess, update displacement and velocity back to elements

    # FIXME: support for variable number of dofs / node
    function to_dict(u, dim, nnodes)
        return Dict(j => [u[dim*(j-1)+k] for k=1:dim] for j=1:nnodes)
    end

    info("Number of nodes = $nnodes, number of time steps = $(length(sol))")

    for (time, x) in zip(sol.t, sol.u)
        full_vec[nz] = x[1:ndofs]
        u = to_dict(full_vec, dim, nnodes)
        full_vec[nz] = x[ndofs+1:end]
        v = to_dict(full_vec, dim, nnodes)
        for problem in get_problems(analysis)
            is_field_problem(problem) || continue
            for elements in get_elements(problem)
                update!(elements, "displacement", time => u)
                update!(elements, "velocity", time => v)
            end
        end
    end

end

# function FEMBase.write_results!(analysis::ImplicitDynamics, writer::Xdmf)
#     # Not implemented yet
# end

export LinearImplicit

end


text:
@testset "heideldiag.jl" begin
    samples = randn(100)

    @testset "results" begin
        @test @inferred(heideldiag(samples)) isa NamedTuple
    end
end


text:
import Core.Compiler

abstract type AbstractThreadPool end

_detect_type(fn, itr) = Core.Compiler.return_type(fn, Tuple{eltype(itr)})

"""
    tforeach(fn, pool, itr)

Mimics `Base.foreach`, but launches the function evaluations onto the provided
pool to assign the tasks.

# Example
```
julia> pool = twith(ThreadPools.LoggedQueuePool(1,2)) do pool
         tforeach(x -> println((x,Threads.threadid())), pool, 1:8)
       end;
(2, 2)
(1, 1)
(3, 2)
(5, 2)
(4, 1)
(6, 2)
(7, 1)
(8, 2)

julia> plot(pool)
```
"""
function tforeach(fn, pool::AbstractThreadPool, itr)
    tmap(fn, pool, itr)
    nothing
end

tforeach(fn, pool::AbstractThreadPool, itr1, itrs...) = tforeach(x -> fn(x...), pool, zip(itr1, itrs...))

"""
    tmap(fn, pool, itr)

Mimics `Base.map`, but launches the function evaluations onto the provided
pool to assign the tasks.

# Example
```
julia> pool = twith(ThreadPools.LoggedQueuePool(1,2)) do pool
         tmap(pool, 1:8) do x
           println((x,Threads.threadid()))
         end
       end;
(2, 2)
(1, 1)
(3, 2)
(4, 1)
(5, 2)
(6, 1)
(7, 2)
(8, 1)

julia> plot(pool)
```
"""
tmap(fn, pool::AbstractThreadPool, itr1, itrs...) = tmap(x -> fn(x...), pool, zip(itr1, itrs...))

"""
    twith(fn, pool) -> pool

Apply the functon `fn` to the provided pool and close the pool.  Returns the
closed pool for any desired analysis or plotting.

# Example
```
julia> twith(ThreadPools.QueuePool(1,2)) do pool
         tforeach(x -> println((x,Threads.threadid())), pool, 1:8)
       end;
(2, 2)
(1, 1)
(3, 2)
(4, 1)
(5, 2)
(6, 1)
(7, 2)
(8, 1)
```
Note in the above example, only two threads were used, as set by the
`QueuePool` setting.
"""
function twith(fn, pool)
    fn(pool)
    close(pool)
    pool
end

"""
    @tthreads pool

Mimic the `Base.Threads.@threads` macro, but uses the provided pool to
assign the tasks.

# Example
```
julia> twith(ThreadPools.QueuePool(1,2)) do pool
         @tthreads pool for x in 1:8
           println((x,Threads.threadid()))
         end
       end;
(2, 2)
(3, 2)
(1, 1)
(4, 2)
(5, 1)
(6, 2)
(8, 2)
(7, 1)
```
"""
macro tthreads(pool, args...)
    na = length(args)
    if na != 1
        throw(ArgumentError("wrong number of arguments in @tthreads"))
    end
    ex = args[1]
    if !isa(ex, Expr)
        throw(ArgumentError("need an expression argument to @tthreads"))
    end
    if ex.head === :for
        if ex.args[1] isa Expr && ex.args[1].head === :(=)
            index = ex.args[1].args[1]
            range = ex.args[1].args[2]
            body = ex.args[2]
            return quote
                tforeach($(esc(pool)), $(esc(range))) do $(esc(index))
                    $(esc(body))
                end
            end
        else
            throw(ArgumentError("nested outer loops are not currently supported by @tthreads"))
        end
    else
        throw(ArgumentError("unrecognized argument to @tthreads"))
    end
end

function Base.finalize(pool::AbstractThreadPool)
    close(pool)
end

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

Closes the pool, shuts down any handlers and finalizes any logging activities.
"""
function Base.close(pool::AbstractThreadPool)
    nothing
end


text:
###############
# API methods #
###############

"""
    ForwardDiff.jacobian(f, x::AbstractArray, cfg::JacobianConfig = JacobianConfig(f, x))

Return `J(f)` evaluated at `x`, assuming `f` is called as `f(x)`.

This method assumes that `isa(f(x), AbstractArray)`.
"""
function jacobian(f, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f, x)) where {T}
    checktag(T, f, x)
    if chunksize(cfg) == length(x)
        return vector_mode_jacobian(f, x, cfg)
    else
        return chunk_mode_jacobian(f, x, cfg)
    end
end

"""
    ForwardDiff.jacobian(f!, y::AbstractArray, x::AbstractArray, cfg::JacobianConfig = JacobianConfig(f!, y, x))

Return `J(f!)` evaluated at `x`,  assuming `f!` is called as `f!(y, x)` where the result is
stored in `y`.
"""
function jacobian(f!, y::AbstractArray, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f!, y, x)) where {T}
    checktag(T, f!, x)
    if chunksize(cfg) == length(x)
        return vector_mode_jacobian(f!, y, x, cfg)
    else
        return chunk_mode_jacobian(f!, y, x, cfg)
    end
end

"""
    ForwardDiff.jacobian!(result::Union{AbstractArray,DiffResult}, f, x::AbstractArray, cfg::JacobianConfig = JacobianConfig(f, x))

Compute `J(f)` evaluated at `x` and store the result(s) in `result`, assuming `f` is called
as `f(x)`.

This method assumes that `isa(f(x), AbstractArray)`.
"""
function jacobian!(result::Union{AbstractArray,DiffResult}, f, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f, x)) where {T}
    checktag(T, f, x)
    if chunksize(cfg) == length(x)
        vector_mode_jacobian!(result, f, x, cfg)
    else
        chunk_mode_jacobian!(result, f, x, cfg)
    end
    return result
end

"""
    ForwardDiff.jacobian!(result::Union{AbstractArray,DiffResult}, f!, y::AbstractArray, x::AbstractArray, cfg::JacobianConfig = JacobianConfig(f!, y, x))

Compute `J(f!)` evaluated at `x` and store the result(s) in `result`, assuming `f!` is
called as `f!(y, x)` where the result is stored in `y`.

This method assumes that `isa(f(x), AbstractArray)`.
"""
function jacobian!(result::Union{AbstractArray,DiffResult}, f!, y::AbstractArray, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f!, y, x)) where {T}
    checktag(T, f!, x)
    if chunksize(cfg) == length(x)
        vector_mode_jacobian!(result, f!, y, x, cfg)
    else
        chunk_mode_jacobian!(result, f!, y, x, cfg)
    end
    return result
end

@inline jacobian(f, x::SArray) = vector_mode_jacobian(f, x)
@inline jacobian(f, x::SArray, cfg::JacobianConfig) = jacobian(f, x)

@inline jacobian!(result::Union{AbstractArray,DiffResult}, f, x::SArray) = vector_mode_jacobian!(result, f, x)
@inline jacobian!(result::Union{AbstractArray,DiffResult}, f, x::SArray, cfg::JacobianConfig) = jacobian!(result, f, x)

#####################
# result extraction #
#####################

@generated function extract_jacobian(::Type{T}, ydual::SArray{SY,VY,DY,M},
                                     x::SArray{SX,VX,DX,N}) where {T,SY,VY,DY,M,SX,VX,DX,N}
    result = Expr(:tuple, [:(partials(T, ydual[$i], $j)) for i in 1:M, j in 1:N]...)
    return quote
        $(Expr(:meta, :inline))
        return SArray{Tuple{M,N}}($result)
    end
end

function extract_jacobian(::Type{T}, ydual::AbstractArray, x::SArray{S,V,D,N}) where {T,S,V,D,N}
    result = similar(ydual, valtype(eltype(ydual)), length(ydual), N)
    return extract_jacobian!(T, result, ydual, N)
end

function extract_jacobian!(::Type{T}, result::AbstractArray, ydual::AbstractArray, n) where {T}
    out_reshaped = reshape(result, length(ydual), n)
    for col in 1:size(out_reshaped, 2), row in 1:size(out_reshaped, 1)
        out_reshaped[row, col] = partials(T, ydual[row], col)
    end
    return result
end

function extract_jacobian!(::Type{T}, result::MutableDiffResult, ydual::AbstractArray, n) where {T}
    extract_jacobian!(T, DiffResults.jacobian(result), ydual, n)
    return result
end

function extract_jacobian_chunk!(::Type{T}, result, ydual, index, chunksize) where {T}
    offset = index - 1
    for i in 1:chunksize
        col = i + offset
        for row in eachindex(ydual)
            result[row, col] = partials(T, ydual[row], i)
        end
    end
    return result
end

reshape_jacobian(result, ydual, xdual) = reshape(result, length(ydual), length(xdual))
reshape_jacobian(result::DiffResult, ydual, xdual) = reshape_jacobian(DiffResults.jacobian(result), ydual, xdual)

###############
# vector mode #
###############

function vector_mode_jacobian(f::F, x, cfg::JacobianConfig{T,V,N}) where {F,T,V,N}
    ydual = vector_mode_dual_eval(f, x, cfg)
    result = similar(ydual, valtype(eltype(ydual)), length(ydual), N)
    extract_jacobian!(T, result, ydual, N)
    extract_value!(T, result, ydual)
    return result
end

function vector_mode_jacobian(f!::F, y, x, cfg::JacobianConfig{T,V,N}) where {F,T,V,N}
    ydual = vector_mode_dual_eval(f!, y, x, cfg)
    map!(d -> value(T,d), y, ydual)
    result = similar(y, length(y), N)
    extract_jacobian!(T, result, ydual, N)
    map!(d -> value(T,d), y, ydual)
    return result
end

function vector_mode_jacobian!(result, f::F, x, cfg::JacobianConfig{T,V,N}) where {F,T,V,N}
    ydual = vector_mode_dual_eval(f, x, cfg)
    extract_jacobian!(T, result, ydual, N)
    extract_value!(T, result, ydual)
    return result
end

function vector_mode_jacobian!(result, f!::F, y, x, cfg::JacobianConfig{T,V,N}) where {F,T,V,N}
    ydual = vector_mode_dual_eval(f!, y, x, cfg)
    map!(d -> value(T,d), y, ydual)
    extract_jacobian!(T, result, ydual, N)
    extract_value!(T, result, y, ydual)
    return result
end

@inline function vector_mode_jacobian(f::F, x::SArray{S,V}) where {F,S,V}
    T = typeof(Tag(f,V))
    return extract_jacobian(T, static_dual_eval(T, f, x), x)
end

@inline function vector_mode_jacobian!(result, f::F, x::SArray{S,V,D,N}) where {F,S,V,D,N}
    T = typeof(Tag(f,V))
    ydual = static_dual_eval(T, f, x)
    result = extract_jacobian!(T, result, ydual, N)
    result = extract_value!(T, result, ydual)
    return result
end

@inline function vector_mode_jacobian!(result::ImmutableDiffResult, f::F, x::SArray{S,V,D,N}) where {F,S,V,D,N}
    T = typeof(Tag(f,V))
    ydual = static_dual_eval(T, f, x)
    result = DiffResults.jacobian!(result, extract_jacobian(T, ydual, x))
    result = DiffResults.value!(d -> value(T,d), result, ydual)
    return result
end

# chunk mode #
#------------#

function jacobian_chunk_mode_expr(work_array_definition::Expr, compute_ydual::Expr,
                                  result_definition::Expr, y_definition::Expr)
    return quote
        @assert length(x) >= N "chunk size cannot be greater than length(x) ($(N) > $(length(x)))"

        # precalculate loop bounds
        xlen = length(x)
        remainder = xlen % N
        lastchunksize = ifelse(remainder == 0, N, remainder)
        lastchunkindex = xlen - lastchunksize + 1
        middlechunks = 2:div(xlen - lastchunksize, N)

        # seed work arrays
        $(work_array_definition)
        seeds = cfg.seeds

        # do first chunk manually to calculate output type
        seed!(xdual, x, 1, seeds)
        $(compute_ydual)
        $(result_definition)
        out_reshaped = reshape_jacobian(result, ydual, xdual)
        extract_jacobian_chunk!(T, out_reshaped, ydual, 1, N)
        seed!(xdual, x, 1)

        # do middle chunks
        for c in middlechunks
            i = ((c - 1) * N + 1)
            seed!(xdual, x, i, seeds)
            $(compute_ydual)
            extract_jacobian_chunk!(T, out_reshaped, ydual, i, N)
            seed!(xdual, x, i)
        end

        # do final chunk
        seed!(xdual, x, lastchunkindex, seeds, lastchunksize)
        $(compute_ydual)
        extract_jacobian_chunk!(T, out_reshaped, ydual, lastchunkindex, lastchunksize)

        $(y_definition)

        return result
    end
end

@eval function chunk_mode_jacobian(f::F, x, cfg::JacobianConfig{T,V,N}) where {F,T,V,N}
    $(jacobian_chunk_mode_expr(quote
                                   xdual = cfg.duals
                                   seed!(xdual, x)
                               end,
                               :(ydual = f(xdual)),
                               :(result = similar(ydual, valtype(eltype(ydual)), length(ydual), xlen)),
                               :()))
end

@eval function chunk_mode_jacobian(f!::F, y, x, cfg::JacobianConfig{T,V,N}) where {F,T,V,N}
    $(jacobian_chunk_mode_expr(quote
                                   ydual, xdual = cfg.duals
                                   seed!(xdual, x)
                               end,
                               :(f!(seed!(ydual, y), xdual)),
                               :(result = similar(y, length(y), xlen)),
                               :(map!(d -> value(T,d), y, ydual))))
end

@eval function chunk_mode_jacobian!(result, f::F, x, cfg::JacobianConfig{T,V,N}) where {F,T,V,N}
    $(jacobian_chunk_mode_expr(quote
                                   xdual = cfg.duals
                                   seed!(xdual, x)
                               end,
                               :(ydual = f(xdual)),
                               :(),
                               :(extract_value!(T, result, ydual))))
end

@eval function chunk_mode_jacobian!(result, f!::F, y, x, cfg::JacobianConfig{T,V,N}) where {F,T,V,N}
    $(jacobian_chunk_mode_expr(quote
                                   ydual, xdual = cfg.duals
                                   seed!(xdual, x)
                               end,
                               :(f!(seed!(ydual, y), xdual)),
                               :(),
                               :(extract_value!(T, result, y, ydual))))
end


text:
immutable PnL
    opened::Union(Date, DateTime)
    closed::Union(Date, DateTime)
    quantity::Int
    commission::Float64
    pnl::Float64
    asset::FinancialAsset
end


text:
# Autogenerated wrapper script for SPIRV_Tools_jll for x86_64-apple-darwin14-cxx11
export spirv_cfg, spirv_reduce, spirv_val, libSPIRV_Tools, spirv_opt, spirv_as, spirv_link, spirv_dis

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "DYLD_FALLBACK_LIBRARY_PATH"
LIBPATH_default = "~/lib:/usr/local/lib:/lib:/usr/lib"

# Relative path to `spirv_cfg`
const spirv_cfg_splitpath = ["bin", "spirv-cfg"]

# This will be filled out by __init__() for all products, as it must be done at runtime
spirv_cfg_path = ""

# spirv_cfg-specific global declaration
function spirv_cfg(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(spirv_cfg_path)
    end
end

# Relative path to `spirv_reduce`
const spirv_reduce_splitpath = ["bin", "spirv-reduce"]

# This will be filled out by __init__() for all products, as it must be done at runtime
spirv_reduce_path = ""

# spirv_reduce-specific global declaration
function spirv_reduce(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(spirv_reduce_path)
    end
end

# Relative path to `spirv_val`
const spirv_val_splitpath = ["bin", "spirv-val"]

# This will be filled out by __init__() for all products, as it must be done at runtime
spirv_val_path = ""

# spirv_val-specific global declaration
function spirv_val(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(spirv_val_path)
    end
end

# Relative path to `libSPIRV_Tools`
const libSPIRV_Tools_splitpath = ["lib", "libSPIRV-Tools-shared.dylib"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libSPIRV_Tools_path = ""

# libSPIRV_Tools-specific global declaration
# This will be filled out by __init__()
libSPIRV_Tools_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libSPIRV_Tools = "@rpath/libSPIRV-Tools-shared.dylib"

# Relative path to `spirv_opt`
const spirv_opt_splitpath = ["bin", "spirv-opt"]

# This will be filled out by __init__() for all products, as it must be done at runtime
spirv_opt_path = ""

# spirv_opt-specific global declaration
function spirv_opt(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(spirv_opt_path)
    end
end

# Relative path to `spirv_as`
const spirv_as_splitpath = ["bin", "spirv-as"]

# This will be filled out by __init__() for all products, as it must be done at runtime
spirv_as_path = ""

# spirv_as-specific global declaration
function spirv_as(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(spirv_as_path)
    end
end

# Relative path to `spirv_link`
const spirv_link_splitpath = ["bin", "spirv-link"]

# This will be filled out by __init__() for all products, as it must be done at runtime
spirv_link_path = ""

# spirv_link-specific global declaration
function spirv_link(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(spirv_link_path)
    end
end

# Relative path to `spirv_dis`
const spirv_dis_splitpath = ["bin", "spirv-dis"]

# This will be filled out by __init__() for all products, as it must be done at runtime
spirv_dis_path = ""

# spirv_dis-specific global declaration
function spirv_dis(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(spirv_dis_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"SPIRV_Tools")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    global spirv_cfg_path = normpath(joinpath(artifact_dir, spirv_cfg_splitpath...))

    push!(PATH_list, dirname(spirv_cfg_path))
    global spirv_reduce_path = normpath(joinpath(artifact_dir, spirv_reduce_splitpath...))

    push!(PATH_list, dirname(spirv_reduce_path))
    global spirv_val_path = normpath(joinpath(artifact_dir, spirv_val_splitpath...))

    push!(PATH_list, dirname(spirv_val_path))
    global libSPIRV_Tools_path = normpath(joinpath(artifact_dir, libSPIRV_Tools_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libSPIRV_Tools_handle = dlopen(libSPIRV_Tools_path)
    push!(LIBPATH_list, dirname(libSPIRV_Tools_path))

    global spirv_opt_path = normpath(joinpath(artifact_dir, spirv_opt_splitpath...))

    push!(PATH_list, dirname(spirv_opt_path))
    global spirv_as_path = normpath(joinpath(artifact_dir, spirv_as_splitpath...))

    push!(PATH_list, dirname(spirv_as_path))
    global spirv_link_path = normpath(joinpath(artifact_dir, spirv_link_splitpath...))

    push!(PATH_list, dirname(spirv_link_path))
    global spirv_dis_path = normpath(joinpath(artifact_dir, spirv_dis_splitpath...))

    push!(PATH_list, dirname(spirv_dis_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
export
    AbstractMAPF,
    MAPF,
    action_type,
    state_type,
    cost_type,
    num_agents,
    num_goals,
    get_env,
    get_starts,
    get_goals,
    get_start,
    get_goal,
    get_start

abstract type AbstractMAPF end

"""
    MAPF

A MAPF is a Multi Agent Path Finding problem. It consists of an environment,
`env`, through which a group of agents may navigate, as well as a list of
start and goal states in that environment. Note that this is the _labeled_
case, where each agent has a specific assigned destination.

Elements:
- env::E - the base environment
- starts::Vector{S} - the vector of initial states
- starts::Vector{G} - the vector of goals
"""
struct MAPF{E,S,G} <: AbstractMAPF # Multi Agent Path Finding Problem
    env     ::E           # Environment Type
    starts  ::Vector{S}   # Vector of initial states
    goals   ::Vector{G}   # Vector of goal states
end
get_env(mapf::MAPF)             = mapf.env
action_type(mapf::MAPF)         = action_type(get_env(mapf))
state_type(mapf::MAPF)          = state_type(get_env(mapf))
cost_type(mapf::MAPF)           = cost_type(get_env(mapf))
num_agents(mapf::MAPF)          = length(mapf.starts)
num_goals(mapf::MAPF)           = length(mapf.goals)
get_starts(mapf::MAPF)          = mapf.starts
get_goals(mapf::MAPF)           = mapf.goals
get_start(mapf::MAPF, i)        = get_starts(mapf)[i]
get_goal(mapf::MAPF, i)         = get_goals(mapf)[i]
get_start(mapf::MAPF, env, i)   = get_start(mapf,i)
base_env_type(mapf::MAPF)       = base_env_type(mapf.env)
base_env_type(env)              = typeof(env)
# TODO implement a check to be sure that no two agents have the same goal
get_initial_cost(mapf::M) where {M<:AbstractMAPF}           = get_initial_cost(mapf.env)
get_infeasible_cost(mapf::M) where {M<:AbstractMAPF}        = get_infeasible_cost(mapf.env)
GraphUtils.get_distance(mapf::AbstractMAPF,args...) = get_distance(mapf.env,args...)


text:
"""
    Tem_F2C(T_degF::Real)
    Tem_C2F(T_degC::Real)
"""
function Tem_F2C(T_degF::Real)
    (T_degF .- 32) ./ (9 / 5)
end

function Tem_C2F(T_degC::Real)
    T_degC .* (9 / 5) .+ 32 # T_degF
end

function heat_index(t::Union{<:Real,Missing}, rh::Union{<:Real,Missing}; missval = -9999)
    if t == missval || rh == missval; return(missing); end
    t = Tem_C2F(t)
    # if (ismissing(t) || ismissing(rh)); return(NaN64); end
    if (t <= 40)
        hi = t
    else
        alpha = 61 + ((t - 68) * 1.2) + (rh * 0.094)
        hi = 0.5 * (alpha + t)
        if (hi > 79)
            hi = -42.379 + 2.04901523 * t + 10.14333127 * rh -
                 0.22475541 * t * rh - 6.83783 * 10^-3 * t^2 -
                 5.481717 * 10^-2 * rh^2 + 1.22874 * 10^-3 * t^2 * rh +
                 8.5282 * 10^-4 * t * rh^2 - 1.99 * 10^-6 * t^2 * rh^2
            if (rh <= 13 && t >= 80 && t <= 112)
                adjustment1 = (13 - rh) / 4
                adjustment2 = sqrt((17 - abs(t - 95)) / 17)
                total_adjustment = adjustment1 * adjustment2
                hi = hi - total_adjustment
            elseif (rh > 85 && t >= 80 && t <= 87)
                adjustment1 = (rh - 85) / 10
                adjustment2 = (87 - t) / 5
                total_adjustment = adjustment1 * adjustment2
                hi = hi + total_adjustment
            end
        end
    end
    Tem_F2C(hi)
end

function heat_index(t::AbstractArray{T}, rh::AbstractArray{T}; missval = -9999) where {T<:Union{Missing,Real}}
    heat_index.(t, rh; missval = missval)
end

"""
    heat_index(t::Union{<:Real,Missing}, rh::Union{<:Real,Missing})
    heat_index(t::AbstractArray{T}, rh::AbstractArray{T}) where {T<:Union{Missing,Real}}
    heat_index(f_tair::AbstractString, f_rh::AbstractString, outfile::AbstractString;
        raw=true, offset = -273.15,
        varname="HI", type = Float32, compress=1,
        overwrite=false)

# Arguments
- `t`: degC
- `rh`: %

# Examples
```julia
heat_index(f_tair, f_rh, outfile; overwrite = false, raw = true, compress = 1)
```

# References
1. https://www.wpc.ncep.noaa.gov/html/heatindex_equationbody.html
2. Kong, D., Gu, X., Li, J., Ren, G., & Liu, J. (2020). Contributions of Global
    Warming and Urbanization to the Intensification of Human‐Perceived Heatwaves
    Over China. Journal of Geophysical Research: Atmospheres, 125(18), 1–16.
    https://doi.org/10.1029/2019JD032175.
"""
function heat_index(f_tair::AbstractString, f_rh::AbstractString, outfile::AbstractString;
    raw=true, offset = -273.15,
    varname="HI", type = Float32, compress=1,
    missval = -9999.0,
    overwrite=false)

    if !isfile(outfile) || overwrite
        isfile(outfile) && rm(outfile)

        println("reading Tair ...")
        arr_tair = nc_read(f_tair; raw=raw) .+ offset

        println("reading RH ...")
        arr_rh = nc_read(f_rh; raw=raw)

        println("calculating HI ...")
        @time arr_HI = heat_index.(arr_tair, arr_rh; missval = missval)

        println("saving ...")
        dims = ncvar_dim(f_tair)
        @time nc_write(arr_HI, outfile, dims; varname=varname, type=type, compress=compress)
    end
end


text:
using PyPlot

function get_X(xs)
    return sum(filter(x->x==2, xs))
end

function random_choice(xs)
    return rand(collect(1:length(xs)))
end

function best_response(xs, i, N, M)
    if xs[i] == 1
        p = get_X(xs)/(N-1)
    else
        p = (get_X(xs) - 1)/(N-1)
    end
    if p * M[1, 2, 1] + (1-p) * M[1, 1, 1] > p * M[2, 2, 1] + (1-p) * M[2, 1, 1]
        return 1
    elseif p * M[1, 2, 1] + (1-p) * M[1, 1, 1] < p * M[2, 2, 1] + (1-p) * M[2, 1, 1]
        return 2
    else
        return rand([1, 2])
    end
end

function sequential_revisions!(xs, N, M, epsilon)
    i = rand(collect(1:length(xs)))
    r = best_response(xs, i, N, M)
    if epsilon > rand()#突然変異
        if 1/2 > rand()
            xs[i] = 3 - r
        else
            xs[i] = r
        end
    else#
        xs[i] = r
    end
end

function simultaneous_revisions!(xs, N, M, epsilon)
    new_xs = copy(xs)
    for i in 1:length(xs)
        r = best_response(xs, i, N, M)
        if epsilon > rand()#突然変異
            if 1/2 > rand()
                new_xs[i] = 3 - r
            else
                new_xs[i] = r
            end
        else#
            new_xs[i] = r
        end
    end
    xs = new_xs
end

function main_sim(xs, M, N, T, epsilon; simul = false)
    Xs = Array(Int, T)
    for t in 1:T
        if simul
            simultaneous_revisions!(xs, N, M, epsilon)
        else
            sequential_revisions!(xs, N, M, epsilon)
        end
        Xs[t] = get_X(xs)
        #println(xs)
    end
    return Xs
end

type KMR
    xs
    M
    N
    T
    epsilon
    Xs
end

KMR(M, N::Int, epsilon) = KMR(ones(Int, N), M, N, 0, epsilon, Array(Int, T))

function simulate!(kmr::KMR, T; init = 0)
    kmr.xs[1:init] = 2
    kmr.Xs = main_sim(kmr.xs, M, N, T, epsilon)
end

function plot_sample_path(kmr)
    plot(kmr.Xs)
end

function plot_empirical_dist(kmr)
    fig, ax = subplots()
    ax[:hist](kmr.Xs)
end

M = Array(Int, (2, 2, 2))
M[1, 1, :] = [4, 4]
M[1, 2, :] = [0, 3]
M[2, 1, :] = [3, 0]
M[2, 2, :] = [2, 2]

N = 20
T = 1000
epsilon = 0.06
"""
loop = 10
for i in 1:loop
    xs = ones(Int, N)#######自由
    Xs = main_sim(xs, M, N, T, epsilon)
    plot(Xs)
end
"""
kmr = KMR(M, N, epsilon)
simulate!(kmr, T)
plot_sample_path(kmr)
plot_empirical_dist(kmr)


text:
include("test_julia.jl")
include("test_hello_world.jl")
include("test_AddCxx.jl")
include("test_StdVector.jl")
include("test_shared_ptr.jl")


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

function Base.copy(::Model)
    return error("Copying nonlinear problems not yet implemented")
end

function Base.show(io::IO, model::Model)
    println(io, "A Nonlinear.Model with:")
    if model.objective !== nothing
    end
    _plural(s, n) = n == 1 ? " 1 $s" : " $n $(s)s"
    println(io, _plural("objective", model.objective !== nothing ? 1 : 0))
    println(io, _plural("parameter", length(model.parameters)))
    println(io, _plural("expression", length(model.expressions)))
    return print(io, _plural("constraint", length(model.constraints)))
end

"""
    set_objective(model::Model, obj)::Nothing

Parse `obj` into a [`Expression`](@ref) and set as the objective
function of `model`.

`obj` must be a type that is supported by [`parse_expression`](@ref).

To remove the objective, pass `nothing`.

## Examples

```julia
model = Model()
x = MOI.VariableIndex(1)
set_objective(model, :(\$x^2 + 1))
set_objective(model, x)
set_objective(model, nothing)
```
"""
function set_objective(model::Model, obj)
    model.objective = parse_expression(model, obj)
    return
end

function set_objective(model::Model, ::Nothing)
    model.objective = nothing
    return
end

"""
    add_expression(model::Model, expr)::ExpressionIndex

Parse `expr` into a [`Expression`](@ref) and add to `model`. Returns an
[`ExpressionIndex`](@ref) that can be interpolated into other input expressions.

`expr` must be a type that is supported by [`parse_expression`](@ref).

## Examples

```julia
model = Model()
x = MOI.VariableIndex(1)
ex = add_expression(model, :(\$x^2 + 1))
set_objective(model, :(sqrt(\$ex)))
```
"""
function add_expression(model::Model, expr)
    push!(model.expressions, parse_expression(model, expr))
    return ExpressionIndex(length(model.expressions))
end

function Base.getindex(model::Model, index::ExpressionIndex)
    return model.expressions[index.value]
end

"""
    add_constraint(
        model::Model,
        func,
        set::Union{
            MOI.GreaterThan{Float64},
            MOI.LessThan{Float64},
            MOI.Interval{Float64},
            MOI.EqualTo{Float64},
        },
    )

Parse `func` and `set` into a [`Constraint`](@ref) and add to `model`. Returns a
[`ConstraintIndex`](@ref) that can be used to delete the constraint or query
solution information.

## Examples

```julia
model = Model()
x = MOI.VariableIndex(1)
c = add_constraint(model, :(\$x^2), MOI.LessThan(1.0))
```
"""
function add_constraint(
    model::Model,
    func,
    set::Union{
        MOI.GreaterThan{Float64},
        MOI.LessThan{Float64},
        MOI.Interval{Float64},
        MOI.EqualTo{Float64},
    },
)
    f = parse_expression(model, func)
    model.last_constraint_index += 1
    index = ConstraintIndex(model.last_constraint_index)
    model.constraints[index] = Constraint(f, set)
    return index
end

"""
    delete(model::Model, c::ConstraintIndex)::Nothing

Delete the constraint index `c` from `model`.

## Examples

```julia
model = Model()
x = MOI.VariableIndex(1)
c = add_constraint(model, :(\$x^2), MOI.LessThan(1.0))
delete(model, c)
```
"""
function delete(model::Model, c::ConstraintIndex)
    delete!(model.constraints, c)
    return
end

function Base.getindex(model::Model, index::ConstraintIndex)
    return model.constraints[index]
end

function MOI.is_valid(model::Model, index::ConstraintIndex)
    return haskey(model.constraints, index)
end

"""
    add_parameter(model::Model, value::Float64)::ParameterIndex

Add a new parameter to `model` with the default value `value`. Returns a
[`ParameterIndex`](@ref) that can be interpolated into other input expressions
and used to modify the value of the parameter.

## Examples

```julia
model = Model()
x = MOI.VariableIndex(1)
p = add_parameter(model, 1.2)
c = add_constraint(model, :(\$x^2 - \$p), MOI.LessThan(0.0))
```
"""
function add_parameter(model::Model, value::Float64)
    push!(model.parameters, value)
    return ParameterIndex(length(model.parameters))
end

function Base.getindex(model::Model, p::ParameterIndex)
    return model.parameters[p.value]
end

function Base.setindex!(model::Model, value::Real, p::ParameterIndex)
    return model.parameters[p.value] = convert(Float64, value)::Float64
end

"""
    register_operator(
        model::Model,
        op::Symbol,
        nargs::Int,
        f::Function,
        [∇f::Function],
        [∇²f::Function],
    )

Register the user-defined operator `op` with `nargs` input arguments in `model`.

## Univariate functions

 * `f(x::T)::T` must be a function that takes a single input argument `x` and
   returns the function evaluated at `x`. If `∇f` and `∇²f` are not provided,
   `f` must support any `Real` input type `T`.
 * `∇f(x::T)::T` is a function that takes a single input argument `x` and
   returns the first derivative of `f` with respect to `x`. If `∇²f` is not
   provided, `∇f` must support any `Real` input type `T`.
 * `∇²f(x::T)::T` is a function that takes a single input argument `x` and
   returns the second derivative of `f` with respect to `x`.

## Multivariate functions

* `f(x::T...)::T` must be a function that takes a `nargs` input arguments `x`
  and returns the function evaluated at `x`. If `∇f` and `∇²f` are not provided,
  `f` must support any `Real` input type `T`.
* `∇f(g::AbstractVector{T}, x::T...)::T` is a function that takes a cache vector
  `g` of length `length(x)`, and fills each element `g[i]` with the partial
  derivative of `f` with respect to `x[i]`.
* `∇²f(H::AbstractMatrix, x::T...)::T` is a function that takes a matrix `H` and
  fills the lower-triangular components `H[i, j]` with the Hessian of `f` with
  respect to `x[i]` and `x[j]` for `i >= j`.

### Notes for multivariate Hessians

 * `H` has `size(H) == (length(x), length(x))`, but you must not access
   elements `H[i, j]` for `i > j`.
 * `H` is dense, but you do not need to fill structural zeros.
"""
function register_operator(model::Model, op::Symbol, nargs::Int, f::Function...)
    return register_operator(model.operators, op, nargs, f...)
end

"""
    evaluate(
        f::AbstractDict,
        model::Model,
        index::ExpressionIndex,
    )

Evaluate the nonlinear expression `index`, where `f[x]` returns the primal value
of decision variable `x::MOI.VariableIndex`.
"""
function evaluate(
    f::AbstractDict,
    model::Model,
    index::ExpressionIndex;
    kwargs...,
)
    return evaluate(f, model, model[index]; kwargs...)
end

"""
    evaluate(
        f::AbstractDict,
        model::Model,
        expr::Expression,
    )

Evaluate the nonlinear expression `expr`, where `f[x]` returns the primal value
of decision variable `x::MOI.VariableIndex`.
"""
function evaluate(
    f::AbstractDict,
    model::Model,
    expr::Expression;
    evaluated_expressions = Dict{Int,Float64}(),
)
    storage = zeros(length(expr.nodes))
    adj = adjacency_matrix(expr.nodes)
    children_arr = SparseArrays.rowvals(adj)
    # An arbitrary limit on the potential input size of a multivariate
    # operation. This will get resized if need-be.
    input_cache = zeros(10)
    for k in length(expr.nodes):-1:1
        node = expr.nodes[k]
        if node.type == NODE_MOI_VARIABLE
            storage[k] = f[MOI.VariableIndex(node.index)]
        elseif node.type == NODE_VALUE
            storage[k] = expr.values[node.index]
        elseif node.type == NODE_SUBEXPRESSION
            if !haskey(evaluated_expressions, node.index)
                evaluated_expressions[node.index] = evaluate(
                    f,
                    model,
                    ExpressionIndex(node.index);
                    evaluated_expressions = evaluated_expressions,
                )
            end
            storage[k] = evaluated_expressions[node.index]
        elseif node.type == NODE_PARAMETER
            storage[k] = model.parameters[node.index]
        elseif node.type == NODE_CALL_MULTIVARIATE
            children_indices = SparseArrays.nzrange(adj, k)
            N = length(children_indices)
            if length(input_cache) < N
                resize!(input_cache, N)
            end
            f_input = view(input_cache, 1:N)
            for (r, i) in enumerate(children_indices)
                f_input[r] = storage[children_arr[i]]
            end
            storage[k] = eval_multivariate_function(
                model.operators,
                model.operators.multivariate_operators[node.index],
                f_input,
            )
        elseif node.type == NODE_CALL_UNIVARIATE
            child_idx = children_arr[adj.colptr[k]]
            storage[k] = eval_univariate_function(
                model.operators,
                model.operators.univariate_operators[node.index],
                storage[child_idx],
            )
        elseif node.type == NODE_COMPARISON
            children_idx = SparseArrays.nzrange(adj, k)
            result = true
            for r in 2:length(children_idx)
                lhs = children_arr[children_idx[r-1]]
                rhs = children_arr[children_idx[r]]
                result &= eval_comparison_function(
                    model.operators,
                    model.operators.comparison_operators[node.index],
                    storage[lhs],
                    storage[rhs],
                )
            end
            storage[k] = result
        else
            @assert node.type == NODE_LOGIC
            children_idx = SparseArrays.nzrange(adj, k)
            lhs = children_arr[children_idx[1]]
            rhs = children_arr[children_idx[2]]
            storage[k] = eval_logic_function(
                model.operators,
                model.operators.logic_operators[node.index],
                storage[lhs] == 1,
                storage[rhs] == 1,
            )
        end
    end
    return storage[1]
end


text:
using ClusterTrees
using StaticArrays

struct Data
    sector::Int
    values::Vector{Int}
end

Data(s=0) = Data(s, Int[])

function sector_center_size(pt, ct, hs)
    hs = hs / 2
    bl = SVector([x > y for (x,y) in zip(pt,ct)]...)
    ct = SVector([b ? y+hs : y-hs for (b,y) in zip(bl,ct)]...)
    sc = sum(b ? 2^(i-1) : 0 for (i,b) in enumerate(bl))
    return sc, ct, hs
end

c = SVector{3,Float64}(1,1,1)
p = SVector{3,Float64}(0,3,-2)

s, c2, h = sector_center_size(p, c, 1.0)
s == 2
c2 == SVector{3,Float64}(0.5,1.5,0.5)
h = 0.5

function lastnonemptychild(tree,node)
    @assert ClusterTrees.PointerBasedTrees.getnode(tree,node).first_child >= 1
    r = 0
    for chd in children(tree, node)
        if !(ClusterTrees.PointerBasedTrees.getnode(tree,chd).first_child < 1)
            r = chd
        end
    end
    return r
end

function lastchild(tree,node)
    @assert ClusterTrees.PointerBasedTrees.getnode(tree,node).first_child >= 1
    r = 0
    for chd in children(tree, node)
        r = chd
    end
    return r
end

struct Router{P,T}
    target_point::P
    smallest_box_size::T
end

function reached(tree, target, meta)
    sector, center, size = meta
    return size <= target.smallest_box_size
end

function directions(tree, target, meta)
    sector, center, size = meta
    return sector_center_size(target.target_point, center, size)
end

function isontherighttrack(tree, node, meta)
    sector, center, size = meta
    ClusterTrees.data(tree, node).sector == sector
end

import ClusterTrees: start, next, done

function ClusterTrees.route!(tree::ClusterTrees.PointerBasedTrees.PointerBasedTree, state, target)

    (parent, prev_fat_par, meta) = state
    reached(tree, target, meta) && return state

    meta = directions(tree, target, meta)
    prev_fat_child = prev_fat_par < 1 ? 0 : lastnonemptychild(tree, prev_fat_par)

    chds = children(tree, parent)
    pos = start(chds)
    while !done(chds, pos)
        child, newpos = next(chds, pos)
        isontherighttrack(tree, child, meta) && return (child, prev_fat_child, meta)
        ClusterTrees.haschildren(tree, child) && (prev_fat_child = child)
        pos = newpos
    end

    sector, center, size = meta
    child = insert!(chds, Data(sector), pos)

    return child, prev_fat_child, meta
end

const N = ClusterTrees.PointerBasedTrees.Node{Data}
tree = ClusterTrees.PointerBasedTrees.PointerBasedTree(
    N[N(Data(), 0, 0, 0, 0)], 1)

function update!(f, tree::ClusterTrees.PointerBasedTrees.PointerBasedTree, i::Int, point, sms::Float64)
    router! = Router(point, sms)
    prev_fat_child = 0
    root_sector, root_center, root_size = 0, SVector{3,Float64}(0,0,0), 1.0
    root_meta = root_sector, root_center, root_size
    root_state = (root(tree), prev_fat_child, root_meta)
    ClusterTrees.update!(f, tree, root_state, i, router!)
end

using DelimitedFiles
Q = readdlm(joinpath(@__DIR__, "points.dlm"), Float64)
points = [SVector{3,Float64}(Q[i,:]) for i in axes(Q,1)]
num_points = length(points)

smallest_box_size = 0.2
f(tree, state, data) = push!(ClusterTrees.data(tree, state[1]).values, data)

for i in 1:num_points
    update!(tree, i, points[i], smallest_box_size) do tree, node, data
        push!(ClusterTrees.data(tree, node).values, data)
    end
end

ns = sum(length(data(tree,nd).values) for nd in ClusterTrees.DepthFirstIterator(tree, tree.root))
@assert ns == num_points

ClusterTrees.print_tree(tree)

num_printed = 0
for node in ClusterTrees.DepthFirstIterator(tree, tree.root)
    global num_printed += 1
end
@assert num_printed == length(tree.nodes)


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BlockPartitioner(sides)
    BlockPartitioner(side₁, side₂, ...)

A method for partitioning spatial objects into blocks of given `sides`.
"""
struct BlockPartitioner{T,N} <: AbstractPartitioner
  sides::SVector{N,T}
end

BlockPartitioner(sides::NTuple{N,T}) where {N,T} =
  BlockPartitioner{T,N}(sides)

BlockPartitioner(sides::Vararg{T,N}) where {N,T} =
  BlockPartitioner(sides)

function partition(object::AbstractSpatialObject{T,N},
                   partitioner::BlockPartitioner{U,N}) where {N,T,U}
  psides = partitioner.sides
  bbox = boundbox(object)

  @assert all(psides .≤ sides(bbox)) "invalid block sides"

  # bounding box properties
  ce = center(bbox)
  lo = lowerleft(bbox)
  up = upperright(bbox)

  # find number of blocks to left and right
  nleft  = @. ceil(Int, (ce - lo) / psides)
  nright = @. ceil(Int, (up - ce) / psides)

  origin  = @. ce - nleft * psides
  nblocks = @. nleft + nright

  subsets = [Vector{Int}() for i in 1:prod(nblocks)]
  neighbors = [Vector{Int}() for i in 1:prod(nblocks)]

  # Cartesian to linear indices
  linear = LinearIndices(Dims(nblocks))

  coords = MVector{N,T}(undef)
  for j in 1:npoints(object)
    coordinates!(coords, object, j)

    # find block coordinates
    c = @. floor(Int, (coords - origin) / psides) + 1
    @inbounds for i in 1:N
      c[i] = clamp(c[i], 1, nblocks[i])
    end
    bcoords = CartesianIndex(Tuple(c))

    # block index
    i = linear[bcoords]

    append!(subsets[i], j)
  end

  # neighboring blocks metadata
  bstart  = CartesianIndex(ntuple(i -> 1, N))
  boffset = CartesianIndex(ntuple(i -> 1, N))
  bfinish = CartesianIndex(Dims(nblocks))
  for (i, bcoords) in enumerate(bstart:bfinish)
    for b in (bcoords - boffset):(bcoords + boffset)
      if all(Tuple(bstart) .≤ Tuple(b) .≤ Tuple(bfinish)) && b ≠ bcoords
        push!(neighbors[i], linear[b])
      end
    end
  end

  # filter out empty blocks
  empty = isempty.(subsets)
  subsets = subsets[.!empty]
  neighbors = neighbors[.!empty]
  for i in findall(empty)
    for n in neighbors
      setdiff!(n, i)
    end
  end

  # save metadata
  metadata = Dict(:neighbors => neighbors)

  SpatialPartition(object, subsets, metadata)
end


text:
abstract type DESPOTRandomSource end

check_consistency(rs::DESPOTRandomSource) = nothing

include("memorizing_rng.jl")

import Random.MT_CACHE_F

mutable struct MemorizingSource{RNG<:AbstractRNG} <: DESPOTRandomSource
    rng::RNG
    memory::Vector{Float64}
    rngs::Matrix{MemorizingRNG{MemorizingSource{RNG}}}
    furthest::Int
    min_reserve::Int
    grow_reserve::Bool
    move_count::Int
    move_warning::Bool
end

function MemorizingSource(K::Int, depth::Int, rng::AbstractRNG; min_reserve=0, grow_reserve=true, move_warning=true)
    RNG = typeof(rng)
    memory = Float64[]
    rngs = Matrix{MemorizingRNG{MemorizingSource{RNG}}}(undef, depth+1, K)
    src = MemorizingSource{RNG}(rng, memory, rngs, 0, min_reserve, grow_reserve, 0, move_warning)
    for i in 1:K
        for j in 1:depth+1
            rngs[j, i] = MemorizingRNG(src.memory, 1, 0, 0, src)
        end
    end
    return src
end

function get_rng(s::MemorizingSource, scenario::Int, depth::Int)
    rng = s.rngs[depth+1, scenario]
    if rng.finish == 0
        rng.start = s.furthest+1
        rng.idx = rng.start - 1
        rng.finish = s.furthest
        reserve(rng, s.min_reserve)
    end
    rng.idx = rng.start - 1
    return rng
end

function Random.seed!(s::MemorizingSource, seed)
    Random.seed!(s.rng, seed)
    resize!(s.memory, 0)
    for i in 1:size(s.rngs, 2)
        for j in 1:size(s.rngs, 1)
            s.rngs[j, i] = MemorizingRNG(s.memory, 1, 0, 0, s)
        end
    end
    s.furthest = 0
    s.move_count = 0
    return s
end

function gen_rand!(r::MemorizingRNG{MemorizingSource{MersenneTwister}}, n::Integer)
    s = r.source

    # make sure that we're on the very end of the source's memory
    if r.finish != s.furthest # we need to move the memory to the end
        len = r.finish - r.start + 1
        if length(s.memory) < s.furthest + len
            resize!(s.memory, s.furthest + len)
        end
        s.memory[s.furthest+1:s.furthest+len] = s.memory[r.start:r.finish]
        r.idx = r.idx - r.start + s.furthest + 1
        r.start = s.furthest + 1
        r.finish = s.furthest + len
        s.furthest += len
        s.move_count += 1
        if s.grow_reserve
            s.min_reserve = max(s.min_reserve, len + n)
        end
    end

    @assert r.finish == s.furthest
    orig = length(s.memory)
    if orig < s.furthest + n
        @assert n <= MT_CACHE_F
        resize!(s.memory, orig+MT_CACHE_F)
        Random.gen_rand(s.rng) # could be faster to use dsfmt_fill_array_close1_open2
        s.memory[orig+1:end] = s.rng.vals
        Random.mt_setempty!(s.rng)
    end
    s.furthest += n
    r.finish += n
    return nothing
end

Random.rng_native_52(s::MemorizingSource) = Random.rng_native_52(s.rng)

function check_consistency(s::MemorizingSource)
    if s.move_count > 0.01*length(s.rngs) && s.move_warning
        @warn("""
             DESPOT's MemorizingSource random number generator had to move the memory locations of the rngs $(s.move_count) times. If this number was large, it may be affecting performance (profiling is the best way to tell).

             To suppress this warning, use MemorizingSource(..., move_warning=false).

             To reduce the number of moves, try using MemorizingSource(..., min_reserve=n) and increase n until the number of moves is low (the final min_reserve was $(s.min_reserve)).
             """)
    end
end


text:
#The code generated need further corrections by you. Execute it line by line and correct the errors.
# Copyright © Vardan Andriasyan, Yauhen Yakimovich, Artur Yakimovich 2015
# MIT license.
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including without
# limitation the rights to use, copy, modify, merge, publish, distribute,
# sublicense, and/or sell copies of the Software, and to permit persons to whom
# the Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

using MatlabCompat

using MatlabCompat.ImageTools

using MatlabCompat.MathTools

using Images,FixedPointNumbers

function graythresh2(img)

  #Check whether the input image is empty
  if isempty(img)
    error("Image is empty");
  end
  #Check whether the input image is gray

  if Images.colorspace(img) != "Gray"
    error("Input Image should be grayscale")
  end

  if isa(raw(img),Array{Uint16,2})
    #Convert image to 8bit and return it's raw values to compute the histogram
    image_array = map(Images.BitShift(Uint8,8),raw(img));
  elseif isa(raw(img),Array{Uint8,2})
    image_array = raw(img);
  else
    warn("Input Image is neither Uint8 or Uint16");
    image_array = raw(img);
  end
  #image_array = raw(img)

  expected_number_of_bins = 2^8;

 (range, counts) = hist(reshape(image_array,length(image_array)),expected_number_of_bins);
  methods(hist)

  #get real number of bins to be used below
  number_of_bins = size(counts)[1];
  #Next the algorithm implementation according to the original paper
  p = counts / sum(counts);
  ω = cumsum(p);
  μ = cumsum(p .* (1:number_of_bins));
  μ_t = μ[end];

  sigma_b_squared = (μ_t * ω - μ).^2 ./ (ω .* (1 - ω));

  # Find the location of the maximum value of sigma_b_squared.
  # The maximum may extend over several bins, so average together the
  # locations.  If maxval is NaN, meaning that sigma_b_squared is all NaN,
  # then return 0.
  maxvalue = maximum(sigma_b_squared);
  if isfinite(maxvalue)
    index = findfirst(sigma_b_squared, maxvalue);

    # Normalize the threshold to the range [0, 1].
    threshold = (index - 1) / (expected_number_of_bins - 1);
  else
    threshold = 0.0;
  end

  return threshold
end

numberOfIterations = 10000;

innerExecutionTimes = float64(zeros(1,numberOfIterations));
outerExecutionTimes = float64(zeros(1,numberOfIterations));
img = MatlabCompat.ImageTools.imread("http://matlabcompat.github.io/img/example.tif");

function simpleJanus()
  threshold = graythresh(img); # compute graysacale threshold using Otsu algorithm
  imgbw =im2bw(img, threshold); # create a binary image based on the graysacel image
  labeledbw = bwlabel(imgbw, 4); # lable each connected object in the image
#   numberOfCells = max(reshape(labeledbw, 1,numel(labeledbw))); # count cells
  #disp(strcat("number of objects:", num2str(numberOfCells)));
end

function forSimpleJanus()

  numberOfIterations = 1000;

  for j=1:numberOfIterations
    simpleJanus();
  end

end

  methods(hist)

for i=1:numberOfIterations
  inner = @elapsed(simpleJanus());
  innerExecutionTimes[i] = inner;
end

# for i=1:numberOfIterations
#   outer =  @elapsed(forSimpleJanus());
#   outerExecutionTimes[i] = outer;
# end

println([median(innerExecutionTimes) std(innerExecutionTimes)])
# print([median(outerExecutionTimes)  std(outerExecutionTimes)])


text:
"""
	readdygraphs(filein)
Load data stored in Dygraphs file format to DataFrame

**Input**
* filein: full name of the input file
* datestring: string defining time format (default: "yyyy/mm/dd HH:MM:SS")

**Output**
* dataframe containing columns + timevector in DateTime format

**Example**
```
data = readdygraphs("../test/input/dygraphs_data.tsf";datestring="yyyymmdd");
```
"""
function readdygraphs(filein::String;datestring::String="yyyy/mm/dd HH:MM:SS")
	# Read header
	channels = Array{SubString{String},1}();
	open(filein,"r") do fid
		# Read channel names
		channels = fid |> readline |> x -> replace(x,r" |\r|\""=>"") |>
						x -> split(x,",") |> x -> x[2:end];
	end
	# read data
	dataall = readdlm(filein,',',skipstart=1,header=false);
	# Convert to DataFrame
	data = DataFrame(datetime=DateTime.(preptimevec.(dataall[:,1]),datestring));
	for (i,v) in enumerate(channels)
		data[Symbol(v)] = convert.(Float64,dataall[:,i+1]);
	end
	return data;
end

"""
	writedygraphs(datain,fileout;decimal)
Write DataFrame to Dygraphs csv format

**Input**
* datain: dataframe containing columns + timevector in DateTime format
* fileout: full name of the output file
* decimal: optional output precision for either all or individual columns. Maximum precision is %.10g!
> output date string is fixed: "yyyy/mm/dd HH:MM:SS"

**Example**
```
datain = DataFrame(temp=[10.,11.,12.,14.],grav=[9.8123,9.9,NaN,9.7],
       datetime=[DateTime(2010,1,1),DateTime(2010,1,2),
           DateTime(2010,1,3),DateTime(2010,1,4)]);
writedygraphs(datain,"../test/output/dygraphs_data.csv",decimal=[1,3]);
```
"""
function writedygraphs(datain::DataFrame,fileout::String;decimal=[4])
	channels,timei = FileTools.findchannels(datain);
	# Round data to output precision
	dataout = round2output(datain,decimal,channels);
	# Write header
	open(fileout,"w") do fid
		@printf(fid,"date");
		for (i,v) in enumerate(channels)
			@printf(fid,",%s",string(v));
		end
		@printf(fid,"\n");
		# Write data
		for i = 1:size(datain,1)
			# Write date + time (do not use Dates.format == very slow)
			@printf(fid,"%04i/%02i/%02i %02i:%02i:%02i",
						Dates.year(datain[timei][i]),
						Dates.month(datain[timei][i]),
						Dates.day(datain[timei][i]),
						Dates.hour(datain[timei][i]),
						Dates.minute(datain[timei][i]),
						Dates.second(datain[timei][i]));
			# add data
			flagval = "NaN";
			for j in channels
				if ismissing(dataout[j][i])
					@printf(fid,",%s",flagval);
				else
					@printf(fid,",%.10g",dataout[j][i]);
				end
			end
			@printf(fid,"\n");
		end
	end
end

"""
	preptimevec(timeval)
Auxilliary function to prepare time vector for DateTime conversion
"""
function preptimevec(timeval)
	if typeof(timeval) == Float64
		return timeval |> x -> round(Int,x) |> string
	elseif typeof(timeval) == Int
		return string(timeval)
	else
		return timeval;
	end
end


text:
using Compat
using FactCheck
using FunctionalCollections

facts("@Persistent constructor macro") do

    context("Persistent Vectors") do
        @fact @Persistent([1, 2, 3]) --> pvec([1, 2, 3])
    end

    context("Persistent Hash Maps") do
        @fact @Persistent(Dict("foo" => 1, "bar" => 2)) --> phmap(("foo", 1), ("bar", 2))
    end

    context("Persistent Set") do
        @fact @Persistent(Set(1, 2, 3, 3)) --> pset(1, 2, 3, 3)
    end

end


text:
"""
    reshape2d(x; dims = ndims(x) - 1)

Reshape `x` into a two-dimensional matrix by joining the first dims dimensions, i.e.
`reshape(x, prod(size(x,i) for i in 1:dims), :)`

`dims=ndims(x)-1` (default) is typically used when turning the output of a 4-D convolution
result into a 2-D input for a fully connected layer.

`dims=1` is typically used when turning the 3-D output of an RNN layer into a 2-D input for
a fully connected layer.

`dims=0` will turn the input into a row vector, `dims=ndims(x)` will turn it into a column
vector.

"""
reshape2d(x; dims=ndims(x)-1)=reshape(x, (dims > 0 ? prod(size(x,i) for i in 1:dims) : 1), :)


text:
"""
    ensure_testcases_installed()

Download `nctiles-testcases` Artifact and decompress files if needed.
"""
function ensure_testcases_installed()
    p=dirname(pathof(NCTiles))
    artifact_toml = joinpath(p, "../Artifacts.toml")
    Pkg.ensure_artifact_installed("NCTILES_TESTCASES",artifact_toml)

    fils=["state_3d_set1.0000000732.data.gz","state_3d_set1.0000001428.data.gz","state_3d_set1.0000002172.data.gz",
    "trsp_3d_set1.0000000732.data.gz","trsp_3d_set1.0000001428.data.gz","trsp_3d_set1.0000002172.data.gz"]
    pth=joinpath(NCTiles.NCTILES_TESTCASES,"diags/")
    for i in 1:length(fils)
        if isfile(joinpath(pth,fils[i]))
            run(`gunzip $(joinpath(pth,fils[i]))`)
        end
    end
end

"""
    getnsteps(d)

Helper function: Determines number of time steps in data d. Input d can be
BinData, NCData, TileData, or an Array.
"""
function getnsteps(d)
    if isa(d,BinData)
        if isa(d.fnames,Array)
            nsteps = length(d.fnames)
        else
            nsteps = 1
        end
    elseif isa(d,NCData)
        ds = Dataset(d.fname)
        dsvar = ds[d.varname]
        if hastimedim(ds,dsvar)
            tdim = dimnames(dsvar)[findtimedim(ds,dsvar)]
            nsteps = length(ds[tdim])
        else
            nsteps = 1
        end
    elseif isa(d,TileData)
        nsteps = getnsteps(d.vals)
    elseif isa(d,Array) && isa(d[1],Array)
        nsteps = length(d)
    else
        nsteps = 1
    end

    return nsteps
end

"""
    checkdims(v0::Array,var::NCvar)

Helper function: Checks that the size of the data about to be written to the file
matches the provided dimensions.
"""
function checkdims(v0,var::NCvar)
    # if v0 is a scalar, 0 dims
    ~isa(v0,Array) ? ndimsv0 = 0 : ndimsv0 = length(size(v0))

    if isa(var.dims,Array)
        if (isa(var.values,Array) && isa(var.values[1],Number))
            dimlist = getfield.(var.dims,:name)
        else
            dimlist = getfield.(var.dims[istimedim.(var.dims).==false],:name)
        end
    else
        dimlist = [string(var.dims)]
    end
    if ndimsv0 != length(dimlist)
        dimlist = join(dimlist,", ")
        error("Size of $(var.name) $(size(v0)) does not match its dimension list: ($dimlist)")
    end

end

"""
    istimedim(d::Union{NCvar,NCDatasets.CFVariable})

Helper function: determines whether d is a time dimension.
"""
function istimedim(d::Union{NCvar,NCDatasets.CFVariable})
    if isa(d,NCvar)
        units = lowercase(d.units)
        longname = get(d.atts,"long_name","")
    else
        units = lowercase(get(d.attrib,"units",""))
        longname = get(d.attrib,"long_name","")
    end
    timeUnits = ["minutes","seconds","hours","days","minute","second","hour","day"]
    return any(occursin.(timeUnits,Ref(units))) || longname=="Time coordinate"
end

"""
    findtimedim(v::Array)

Helper function: finds which dimension is a time dimension, if any. Can be Array of
    NCvars or NCDatasets.CFVariables.
"""
function findtimedim(dims::Array)
    return findall(istimedim.(dims))[1]
end

findtimedim(v::NCvar) = findtimedim(v.dims)

function findtimedim(ds::NCDatasets.Dataset,v::NCDatasets.CFVariable)
    idx = [haskey(ds,dim) for dim in dimnames(v)]  # tcb dim isn't always listed as a variable in NetCDF files
    findall(idx)[findtimedim([ds[d] for d in dimnames(v)[idx]])]
end

"""
    hastimedim::Array)

Helper function: determines whether an array of dimensions has a time dimension. Can
    be Array of NCvars or NCDatasets.CFVariables.
"""
function hastimedim(dims::Array)
    return any(istimedim.(dims))
end

"""
    hastimedim(v::NCvar)

Helper function: determines whether a variable has a time dimension.
"""
function hastimedim(v::NCvar)
    ncvardim = isa.(v.dims,NCvar)
    return any(ncvardim) && hastimedim(v.dims[ncvardim])
end

function hastimedim(ds::NCDatasets.Dataset,v::NCDatasets.CFVariable)
    vardims = dimnames(v)[[haskey(ds,dim) for dim in dimnames(v)]] # tcb dim isn't always listed as a variable in NetCDF files
    hastimedim([ds[d] for d in vardims])
end

"""
    getdims(myflds::Dict)

Helper function: retrieves unique dimensions in dictionary of fields.
"""
function getdims(myflds::Dict)
    dims = unique(vcat([myflds[v].dims for v in keys(myflds)]...))
    dims = filter( d -> isa(d,NCvar),dims)
end

"""
    hastiledata(myfld::NCvar)

Helper function: returns whether given field or dict of fields has any tiled data.
"""
function hastiledata(myfld::NCvar)
    return isa(myfld.values,TileData)
end

hastiledata(myflds::Dict) = any([hastiledata(myflds[f]) for f in keys(myflds)])


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

using Test

using QuantumCircuits
using QuantumCircuits.Circuit: getCode
using QuantumCircuits.QML

# Test the graph u3 optimization
expqc = QCircuit(1)
expqc.u3(0, ParameterT(1.0), ParameterT(2.0), ParameterT(3.0))

qc = QCircuit(1)
qc.u3(0, ParameterT(1.0), ParameterT(2.0), ParameterT(3.0))
qc.u3(0, ParameterT(4.0), ParameterT(5.0), ParameterT(6.0))
@test qc == expqc
@test length(getCode(qc)) == 1

# with random parameters
qc = QCircuit(1)
qc.u3(0)
qc.u3(0)
@test length(getparameters(qc)) == 3
@test length(getCode(qc)) == 1

# No optimization in case wthout parameters
qc = QCircuit(1)
qc.u3(0, 1.0, 2.0, 3.0)
qc.u3(0, 4.0, 5.0, 6.0)
@test length(getparameters(qc)) == 0
@test length(getCode(qc)) == 2


text:
module CreateTableFoos

import SearchLight.Migrations: create_table, column, columns, primary_key, add_index, drop_table, add_indices

function up()
  create_table(:foos) do
    [
      primary_key()
      column(:column_name, :column_type)
      columns([
        :column_name => :column_type
      ])
    ]
  end

  add_index(:foos, :column_name)
  add_indices(foos, :column_name_1, :column_name_2)
end

function down()
  drop_table(:foos)
end

end


text:
function mle_for_θ(y::AbstractVector, μ::AbstractVector, wts::AbstractVector;
                   maxIter=30, convTol=1.e-6)
    function first_derivative(θ::Real)
        tmp(yi, μi) = (yi+θ)/(μi+θ) + log(μi+θ) - 1 - log(θ) - digamma(θ+yi) + digamma(θ)
        unit_weights ? sum(tmp(yi, μi) for (yi, μi) in zip(y, μ)) :
                       sum(wti * tmp(yi, μi) for (wti, yi, μi) in zip(wts, y, μ))
    end
    function second_derivative(θ::Real)
        tmp(yi, μi) = -(yi+θ)/(μi+θ)^2 + 2/(μi+θ) - 1/θ - trigamma(θ+yi) + trigamma(θ)
        unit_weights ? sum(tmp(yi, μi) for (yi, μi) in zip(y, μ)) :
                       sum(wti * tmp(yi, μi) for (wti, yi, μi) in zip(wts, y, μ))
    end

    unit_weights = length(wts) == 0
    if unit_weights
        n = length(y)
        θ = n / sum((yi/μi - 1)^2 for (yi, μi) in zip(y, μ))
    else
        n = sum(wts)
        θ = n / sum(wti * (yi/μi - 1)^2 for (wti, yi, μi) in zip(wts, y, μ))
    end
    δ, converged = one(θ), false

    for t = 1:maxIter
        θ = abs(θ)
        δ = first_derivative(θ) / second_derivative(θ)
        if abs(δ) <= convTol
            converged = true
            break
        end
        θ = θ - δ
    end
    converged || throw(ConvergenceException(maxIter))
    θ
end

function negbin(F, D, args...;
                initialθ::Real=Inf, maxIter::Integer=30, convTol::Real=1.e-6,
                verbose::Bool=false, kwargs...)
    maxIter >= 1 || throw(ArgumentError("maxIter must be positive"))
    convTol > 0  || throw(ArgumentError("convTol must be positive"))
    initialθ > 0 || throw(ArgumentError("initialθ must be positive"))

    # fit a Poisson regression model if the user does not specify an initial θ
    if isinf(initialθ)
        regmodel = glm(F, D, Poisson(), args...;
                       maxIter=maxIter, convTol=convTol, verbose=verbose, kwargs...)
    else
        regmodel = glm(F, D, NegativeBinomial(initialθ), args...;
                       maxIter=maxIter, convTol=convTol, verbose=verbose, kwargs...)
    end

    μ = regmodel.model.rr.mu
    y = regmodel.model.rr.y
    wts = regmodel.model.rr.wts
    lw, ly = length(wts), length(y)
    if lw != ly && lw != 0
        throw(ArgumentError("length of wts must be either $ly or 0 but was $lw"))
    end

    θ = mle_for_θ(y, μ, wts)
    d = sqrt(2 * max(1, deviance(regmodel)))
    δ = one(θ)
    ll = loglikelihood(regmodel)
    ll0 = ll + 2 * d

    converged = false
    for i = 1:maxIter
        if abs(ll0 - ll)/d + abs(δ) <= convTol
            converged = true
            break
        end
        verbose && println("[ Alternating iteration ", i, ", θ = ", θ, " ]")
        regmodel = glm(F, D, NegativeBinomial(θ), args...;
                       maxIter=maxIter, convTol=convTol, verbose=verbose, kwargs...)
        μ = regmodel.model.rr.mu
        prevθ = θ
        θ = mle_for_θ(y, μ, wts; maxIter=maxIter, convTol=convTol)
        δ = prevθ - θ
        ll0 = ll
        ll = loglikelihood(regmodel)
    end
    converged || throw(ConvergenceException(maxIter))
    regmodel
end


text:
using MCVI
using Test
import MCVI: init_lower_action, lower_bound, upper_bound
using POMDPs
using Random
using POMDPModelTools
using POMDPSimulators
using POMDPModels # for LightDark1d

# to get tests to pass on POMDPs 0.7.3 - can be removed once support is dropped for 0.7.3
POMDPs.generate_o(m::LightDark1D, sp, rng::AbstractRNG) = rand(rng, observation(m, sp))

# Bounds
mutable struct LightDark1DLowerBound
    rng::AbstractRNG
end

mutable struct LightDark1DUpperBound
    rng::AbstractRNG
end

function lower_bound(lb::LightDark1DLowerBound, p::LightDark1D, s::LightDark1DState)
    _, _, r = gen(DDNOut(:sp,:o,:r), p, s, init_lower_action(p), lb.rng)
    return r * discount(p)
end

function upper_bound(ub::LightDark1DUpperBound, p::LightDark1D, s::LightDark1DState)
    steps = abs(s.y)/p.step_size + 1
    return p.correct_r*(discount(p)^steps)
end

function init_lower_action(p::LightDark1D)
    return 0 # Worst? This depends on the initial state? XXX
end

include("test_policy.jl")

include("test_updater.jl")

include("test_belief.jl")

@test test_dummy_graph()
@test test_dummy_graph2()

@test test_backup()

@test test_belief()

include("test_solve.jl")
@test test_solve()

include("test_simulation.jl")
test_simulation()

include("test_requirements.jl")
test_requirements()

# @test test_dummy_heuristics()


text:
using MOBC, vOptGeneric, CPLEX, PyPlot, Suppressor
using Base.Test

function benchmark( ; fig1 = 1, fig2 = 2, args...)
	res_stidsen=Dict{Int, Vector{Float64}}() ; nb_nodes_stidsen=Dict{Int, Vector{Int}}()
	res_parragh=Dict{Int, Vector{Float64}}() ; nb_nodes_parragh=Dict{Int, Vector{Int}}()
	res_eps=Dict{Int, Vector{Float64}}()
	for filename in ("gap1.txt",)# "gap2.txt")
		filepath = joinpath(Pkg.dir("MOBC"), "test", "Instances_gap", "gap_orlib", filename)
		instances = parseGAP_orlib(filepath)
		n = instances[1].m.numCols
		res_stidsen[n]=Float64[] ; res_parragh[n]=Float64[] ; res_eps[n]=Float64[]
		nb_nodes_stidsen[n]=Int[] ; nb_nodes_parragh[n]=Int[]
		for inst in instances
			m = inst.m
			global m_current = m
			valBC, tBC, _ = @timed solve_stidsen(m, Inf ; args...);
			sleep(0.5)
			valPAR, tPAR, _ = @timed solve_parragh(m, Inf ; args...);
			sleep(0.5)
			status,tEPS,_ = @suppress @timed solve(m, method=:epsilon, round_results=true);
			if !(length(getY_N(m)) == length(unique(valBC.YN)) == length(unique(valPAR.YN)))
				print(m)
				@show valBC.YN
				@show valPAR.YN
				@show getY_N(m)
				figure(3)
				clf()
				plot(map(first, valPAR.YN), map(last,valPAR.YN), "b>", markersize="5", label="parragh")
				plot(map(first, valBC.YN), map(last, valBC.YN), "r<", markersize="5", label="stidsen")
				plot(map(first, getY_N(m)), map(last, getY_N(m)), "go", markersize="2", label="epsilon")
				legend()
				error()
			end

			push!(res_stidsen[n], tBC) ; push!(res_parragh[n], tPAR) ; push!(res_eps[n], tEPS)
			push!(nb_nodes_parragh[n], valPAR.nodes) ; push!(nb_nodes_stidsen[n], valBC.nodes)

			figure(fig1) ; clf()
			xlabel("n") ; ylabel("t(s)") ; title("time per instance")
			plot([],[],"gx",markersize=3,label="stidsen") ; plot([],[],"bx",markersize=3,label="parragh") ; plot([],[],"k.",markersize=3,label="epsilon")
			for k in keys(res_stidsen)
				plot(fill(k, length(res_stidsen[k])), res_stidsen[k], "gx", markersize="3")
				plot(fill(k, length(res_stidsen[k])), res_parragh[k], "bx", markersize="3")
				plot(fill(k, length(res_stidsen[k])), res_eps[k], "k.", markersize="3")
			end
			legend()

			figure(fig2) ; clf()
			xlabel("n") ; ylabel("#nodes") ; title("average time and #nodes")
			k = sort(collect(keys(res_stidsen)))
			plot(k, [mean(nb_nodes_stidsen[a]) for a in k], "gs", markersize="3") ; plot(k, [mean(nb_nodes_stidsen[a]) for a in k], "g--", linewidth=1)
			plot(k, [mean(nb_nodes_parragh[a]) for a in k], "bs", markersize="3") ; plot(k, [mean(nb_nodes_parragh[a]) for a in k], "b--", linewidth=1)
			twinx() ; ylabel("t(s)")
			plot([],[],"gs",markersize="3",label="#nodes stidsen") #labels dont work after twinx(), had to put a dummy plot
			plot([],[],"bs",markersize="3",label="#nodes parragh")
			plot(k, [mean(res_stidsen[a]) for a in k], "gx", markersize="3", label="tmoy stidsen") ; plot(k, [mean(res_stidsen[a]) for a in k], "g-", linewidth=1)
			plot(k, [mean(res_parragh[a]) for a in k], "bx", markersize="3", label="tmoy parragh") ; plot(k, [mean(res_parragh[a]) for a in k], "b-", linewidth=1)
			plot(k, [mean(res_eps[a])for a in k], "b.", markersize="3", label="tmoy epsilon")
			legend()
		end
		println("\n####################################\nn = $n")
		println("tmoy eps : $(mean(res_eps[n]))") ; println("tmoy stidsen : $(mean(res_stidsen[n]))") ; println("tmoy parragh : $(mean(res_parragh[n]))")
		println("moy : $(round(mean(nb_nodes_stidsen[n]),1)) noeuds par BC Stidsen") ; println("moy : $(round(mean(nb_nodes_parragh[n]),1)) noeuds par BC Parragh")
	end
	res_stidsen, res_parragh, res_eps, nb_nodes_stidsen, nb_nodes_parragh
end

srand(0)
benchmark(use_nsga=true, global_branch=false, docovercuts=true, lift_covers=true)


text:
@testset "readers_writer_lock" begin
    rwl = ReaderCountRWLock()

    @test islocked(rwl) == false
    @test is_read_locked(rwl) == false

    read_lock(rwl)
    @test is_read_locked(rwl) == true
    @test islocked(rwl) == false

    read_unlock(rwl)
    @test is_read_locked(rwl) == false
    @test islocked(rwl) == false

    lock(rwl)
    @test is_read_locked(rwl) == false
    @test islocked(rwl) == true

    unlock(rwl)
    @test is_read_locked(rwl) == false
    @test islocked(rwl) == false

    @test_throws ErrorException read_unlock(rwl)
end


text:
using BenchmarkTools
using Plots

const SUITE = BenchmarkGroup()
julia_cmd = split(get(ENV, "TESTCMD", Base.JLOptions().julia_bin))

SUITE["load_plot_display"] = @benchmarkable run(`$julia_cmd --startup-file=no --project -e 'using Plots; display(plot(1:0.1:10, sin.(1:0.1:10)))'`)
SUITE["load"] = @benchmarkable run(`$julia_cmd --startup-file=no --project -e 'using Plots'`)
SUITE["plot"] = @benchmarkable p = plot(1:0.1:10, sin.(1:0.1:10))
SUITE["display"] = @benchmarkable display(p) setup=(p = plot(1:0.1:10, sin.(1:0.1:10)))


text:
using FTCTests
using Test

@test FTCTests.TRAJ_DATA_NAME == "traj.jld2"


text:
# erste Schuss
function findminnorm(list, ego)
    nrpoints = size(list)[2]
    distances = Array{Float64}(nrpoints)
    for i in 1:nrpoints
        distances[i] = norm(list[:,i]-ego)
    end
    findmin(distances)
end

# vektorisiert
function fnp(list, ego)
    list = list.-ego
    list[1,:] = list[1,:].^2
    list[2,:] = list[2,:].^2
    distance = list[1,:]+list[2,:]
    findmin(distance) #returns squared distance
end

function pointsvecnorm(list, ego)
    list2 = map(x-> norm(x-ego), list)
    findmin(list2)
end

function pointsvecnorm_opt(list, ego)
    nrpoints = length(list)
    distances = Array{Float64}(nrpoints)
    for i in 1:nrpoints
        distances[i] = norm(list[i]-ego)
    end
    findmin(distances)
end

max_stars_repo_name: Affos/Bamberg

max_issues_repo_name: Affos/Bamberg

max_forks_repo_name: Affos/Bamberg


text:
addprocs(4)
using ClusterDicts
using Base.Test

function test_d(d, pids, name, ktype, vtype)
    common_tests(d, pids, name, ktype, vtype)

    if isa(d, GlobalDict)
        # test it is the same everywhere
        (k, v) = getkv(100, ktype, vtype)
        d[k] = v
        (_, v) = getkv(101, ktype, vtype)
        d[k] = v
        @sync begin
            results=Bool[]
            for p in d.pids
                @async push!(results, v == remotecall_fetch((n, k)->ClusterDicts.directory[n][k], p, d.name, k))
            end
        end
        @test all(results)

        (k, _) = getkv(102, ktype, vtype)
        v = vtype in (ASCIIString, Any) ? ValueF(()->string(myid())) : ValueF(()->convert(vtype, myid()))
        d[k] = v

        @sync begin
            results=Bool[]
            for p in d.pids
                if vtype in (ASCIIString, Any)
                    @async push!(results, remotecall_fetch((D, K)-> myid() == convert(Int, parse(Float64, D[K])), p, d, k))
                else
                    @async push!(results, remotecall_fetch((D, K)-> myid() == convert(Int, D[K]), p, d, k))
                end
            end
        end
        @test all(results)

    else
        # test it is distributed properly
        for i in 1001:1100
            (k, v) = getkv(i, ktype, vtype)
            d[k] = v
        end

        @sync begin
            results=[]
            for p in d.pids
                @async push!(results, remotecall_fetch(n->length(ClusterDicts.directory[n]), p, d.name))
            end
        end
        @test sum(results) == length(d)

        results = Int[]
        for i in 1200:1299
            (k, _) = getkv(i, ktype, vtype)
            v = vtype in (ASCIIString, Any) ? ValueF(()->string(myid())) : ValueF(()->convert(vtype, myid()))
            d[k] = v

            if vtype in (ASCIIString, Any)
                push!(results, convert(Int, parse(Float64, d[k])))
            else
                push!(results, convert(Int, d[k]))
            end
        end

        (r, counts) = hist(results)
        for c in counts
           @test (c > 1 && c < 100)
        end
        @test convert(Int, maximum(r)) == maximum(d.pids)
    end
end

function getkv(i, ktype, vtype)
    k = ktype in (ASCIIString, Any) ? "$i" : convert(ktype, i)
    v = vtype in (ASCIIString, Any) ? "$i" : convert(vtype, i)
    (k, v)
end

function common_tests(d, pids, name, ktype, vtype)
    println("Testing : ", d)
    @test isempty(d) == true
    (k, v) = getkv(0, ktype, vtype)
    d[k] = v
    @test isempty(d) == false
    @test d[k] == v
    @test length(d) == 1

    (k, v) = getkv(1, ktype, vtype)
    d[k] = v
    @test length(d) == 2
    @test d[k] == v

    (k, v) = getkv(2, ktype, vtype)

    @test_throws Exception d[k]
    @test get(d, k, 1234) == 1234
    @test get!(d, k, v) == v
    @test get(d, k, 1234) == v

    newv = vtype in (ASCIIString, Any) ? "3" : convert(vtype, 3)
    d[k] = newv
    @test d[k] == newv

    if ktype != Any && vtype != Any
        @test_throws Exception gd[:TYPEERR] = v
        @test_throws Exception gd[k] = :TYPEERR
        @test_throws Exception gd[:TYPEERR] = :TYPEERR
    end

end

for pids in Any[PidsAll(), PidsWorkers(), [2,3], :default]
    for KV in Any[(Int, Int), (Int, Float64), (ASCIIString, ASCIIString), :default]
        for name in Any["foobar", :default]
            kwargs=[]
            name != :default && push!(kwargs, (:name, name))
            if KV != :default
                push!(kwargs, (:ktype, KV[1]))
                push!(kwargs, (:vtype, KV[2]))
            else
                KV = (Any, Any)
            end

            for dtype in [GlobalDict, DistributedDict]
                if pids == :default
                    d = dtype(; kwargs...)
                else
                    d = dtype(pids; kwargs...)
                end
                test_d(d, pids, name, KV[1], KV[2])
                delete!(d)
            end
        end
    end
end

println("Testing get/set from both participating and non-participating workers")
println()
for dtype in [GlobalDict, DistributedDict]
    d = dtype([2,3]) # dict on subset of workers

    i = 1
    # test get/set from both participating and non-participating workers
    for p in procs()
        i = i+1; k = i; v = string(i)
        d[k] = v
        @test v == remotecall_fetch((D, K)->D[K], p, d, k)

        i = i+1; k = i; v = string(i)
        remotecall_fetch((D, K, V)->(D[K]=V; nothing), p, d, k, v)
        @test d[k] == v
    end

    delete!(d)
end

# ensure that the directory is empty on all nodes
println("Test for any leaks")
println()
results=Bool[]
for p in procs()
    push!(results, remotecall_fetch(()->isempty(ClusterDicts.directory), p))
end

@test all(results)


text:
@testset "235.lowest-common-ancestor-of-a-binary-search-tree.jl" begin
    root = TreeNode{Int}([6, 2, 8, 0, 4, 7, 9, nothing, nothing, 3, 5])
    @test lowest_common_ancestor_235(root, root.left, root.right) == root
    @test lowest_common_ancestor_235(root, root.left, root.left.right) == root.left
    root = TreeNode{Int}([2, 1])
    @test lowest_common_ancestor_235(root, root, root.left) == root
end


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

mutable struct InstanceViewStatus <: SwaggerModel
    code::Any # spec type: Union{ Nothing, String } # spec name: code
    level::Any # spec type: Union{ Nothing, String } # spec name: level
    displayStatus::Any # spec type: Union{ Nothing, String } # spec name: displayStatus
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    time::Any # spec type: Union{ Nothing, DateTime } # spec name: time

    function InstanceViewStatus(;code=nothing, level=nothing, displayStatus=nothing, message=nothing, time=nothing)
        o = new()
        validate_property(InstanceViewStatus, Symbol("code"), code)
        setfield!(o, Symbol("code"), code)
        validate_property(InstanceViewStatus, Symbol("level"), level)
        setfield!(o, Symbol("level"), level)
        validate_property(InstanceViewStatus, Symbol("displayStatus"), displayStatus)
        setfield!(o, Symbol("displayStatus"), displayStatus)
        validate_property(InstanceViewStatus, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(InstanceViewStatus, Symbol("time"), time)
        setfield!(o, Symbol("time"), time)
        o
    end
end # type InstanceViewStatus

const _property_map_InstanceViewStatus = Dict{Symbol,Symbol}(Symbol("code")=>Symbol("code"), Symbol("level")=>Symbol("level"), Symbol("displayStatus")=>Symbol("displayStatus"), Symbol("message")=>Symbol("message"), Symbol("time")=>Symbol("time"))
const _property_types_InstanceViewStatus = Dict{Symbol,String}(Symbol("code")=>"String", Symbol("level")=>"String", Symbol("displayStatus")=>"String", Symbol("message")=>"String", Symbol("time")=>"DateTime")
Base.propertynames(::Type{ InstanceViewStatus }) = collect(keys(_property_map_InstanceViewStatus))
Swagger.property_type(::Type{ InstanceViewStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_InstanceViewStatus[name]))}
Swagger.field_name(::Type{ InstanceViewStatus }, property_name::Symbol) =  _property_map_InstanceViewStatus[property_name]

const _allowed_InstanceViewStatus_level = ["Info", "Warning", "Error"]

function check_required(o::InstanceViewStatus)
    true
end

function validate_property(::Type{ InstanceViewStatus }, name::Symbol, val)
    if name === Symbol("level")
        Swagger.validate_param(name, "InstanceViewStatus", :enum, val, _allowed_InstanceViewStatus_level)
    end
end