text:
# definitions related to C interface

import Core.Intrinsics.cglobal

# constants to match JL_RTLD_* in src/julia.h
const RTLD_LOCAL     = 0x00000000
const RTLD_GLOBAL    = 0x00000001
const RTLD_LAZY      = 0x00000002
const RTLD_NOW       = 0x00000004
const RTLD_NODELETE  = 0x00000008
const RTLD_NOLOAD    = 0x00000010
const RTLD_DEEPBIND  = 0x00000020
const RTLD_FIRST     = 0x00000040

dlsym(hnd, s::Union(Symbol,String)) = ccall(:jl_dlsym, Ptr{Void}, (Ptr{Void}, Ptr{Uint8}), hnd, s)
dlsym_e(hnd, s::Union(Symbol,String)) = ccall(:jl_dlsym_e, Ptr{Void}, (Ptr{Void}, Ptr{Uint8}), hnd, s)
dlopen(s::String, flags::Integer) = ccall(:jl_load_dynamic_library, Ptr{Void}, (Ptr{Uint8},Uint32), s, flags)
dlopen_e(s::String, flags::Integer) = ccall(:jl_load_dynamic_library_e, Ptr{Void}, (Ptr{Uint8},Uint32), s, flags)
dlopen(s::String) = dlopen(s, RTLD_LAZY | RTLD_DEEPBIND)
dlopen_e(s::String) = dlopen_e(s, RTLD_LAZY | RTLD_DEEPBIND)
dlclose(p::Ptr) = if p!=C_NULL; ccall(:uv_dlclose,Void,(Ptr{Void},),p); end

cfunction(f::Function, r, a) =
    ccall(:jl_function_ptr, Ptr{Void}, (Any, Any, Any), f, r, a)

if ccall(:jl_is_char_signed, Any, ())
    typealias Cchar Int8
else
    typealias Cchar Uint8
end
typealias Cuchar Uint8
typealias Cshort Int16
typealias Cushort Uint16
typealias Cint Int32
typealias Cuint Uint32
if OS_NAME === :Windows
    typealias Clong Int32
    typealias Culong Uint32
    typealias Cwchar_t Uint16
else
    typealias Clong Int
    typealias Culong Uint
    typealias Cwchar_t Int32
end
typealias Cptrdiff_t Int
typealias Csize_t Uint
typealias Cssize_t Int
typealias Clonglong Int64
typealias Culonglong Uint64
typealias Cfloat Float32
typealias Cdouble Float64
#typealias Ccomplex_float Complex64
#typealias Ccomplex_double Complex128

const sizeof_off_t = ccall(:jl_sizeof_off_t, Cint, ())

if sizeof_off_t === 4
    typealias FileOffset Int32
else
    typealias FileOffset Int64
end

typealias Coff_t FileOffset

# deferring (or un-deferring) ctrl-c handler for external C code that
# is not interrupt safe (see also issue #2622).  The sigatomic_begin/end
# functions should always be called in matched pairs, ideally via:
#            disable_sigint() do .. end
# reennable_sigint is provided so that immediate ctrl-c handling is
# re-enabled within a sigatomic region, e.g. inside a Julia callback function
# within a long-running C routine.
sigatomic_begin() = ccall(:jl_sigatomic_begin, Void, ())
sigatomic_end() = ccall(:jl_sigatomic_end, Void, ())
disable_sigint(f::Function) = try sigatomic_begin(); f(); finally sigatomic_end(); end
reenable_sigint(f::Function) = try sigatomic_end(); f(); finally sigatomic_begin(); end

# flush C stdio output from external libraries
flush_cstdio() = ccall(:jl_flush_cstdio, Void, ())

function find_library{T<:ByteString, S<:ByteString}(libnames::Array{T,1}, extrapaths::Array{S,1}=ASCIIString[])
    for lib in libnames
        for path in extrapaths
            l = joinpath(path, lib)
            p = dlopen_e(l, RTLD_LAZY)
            if p != C_NULL
                dlclose(p)
                return l
            end
        end
        p = dlopen_e(lib, RTLD_LAZY)
        if p != C_NULL
            dlclose(p)
            return lib
        end
    end
    return ""
end

function ccallable(f::Callable, rt::Type, argt::(Type...), name::Union(String,Symbol)=string(f))
    ccall(:jl_extern_c, Void, (Any, Any, Any, Ptr{Uint8}), f, rt, argt, name)
end

function ccallable(f::Callable, argt::(Type...), name::Union(String,Symbol)=string(f))
    ccall(:jl_extern_c, Void, (Any, Ptr{Void}, Any, Ptr{Uint8}), f, C_NULL, argt, name)
end

macro ccallable(def)
    if isa(def,Expr) && (def.head === :(=) || def.head === :function)
        sig = def.args[1]
        if sig.head === :call
            name = sig.args[1]
            at = map(sig.args[2:end]) do a
                if isa(a,Expr) && a.head === :(::)
                    a.args[2]
                else
                    :Any
                end
            end
            return quote
                $(esc(def))
                ccallable($(esc(name)), $(Expr(:tuple, map(esc, at)...)))
            end
        end
    end
    error("expected method definition in @ccallable")
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
function champs_2!(champu   :: Array{ComplexF64,2},
                   champv   :: Array{ComplexF64,2},
                   m        :: MicMac,
                   t        :: Float64,
                   fft_ubar :: Vector{ComplexF64},
                   fft_vbar :: Vector{ComplexF64},
                   fft_ug   :: Vector{ComplexF64},
                   fft_vg   :: Vector{ComplexF64})

    ftau!(champu, champv, m, t, fft_ubar, fft_vbar)

    fft!(champu, 1)
    fft!(champv, 1)

    champu[1, :] .= 0.0
    champv[1, :] .= 0.0

    dtauh1u = copy(champu)
    dtauh1v = copy(champv)

    ifft!(dtauh1u, 1)
    ifft!(dtauh1v, 1)

    champu ./= (1im * m.ktau)
    champv ./= (1im * m.ktau)

    ifft!(champu, 1)
    ifft!(champv, 1)

    champu .*= m.epsilon
    champv .*= m.epsilon

    champu .+= transpose(fft_ubar)
    champv .+= transpose(fft_vbar)

    ffu    = copy(champu)
    ffv    = copy(champv)

    ffu   .+= transpose(fft_ug)
    ffv   .+= transpose(fft_vg)

    ftau!(m, t, ffu, ffv)

    ftau!(m, t, champu, champv)

    fft!(champu, 1)
    fft!(champv, 1)

    champubar = champu[1, :] / m.ntau
    champvbar = champv[1, :] / m.ntau

    duftau!(champu, champv, m, t, fft_ubar, fft_vbar, champubar, champvbar)

    dtftau!(champu, champv, m, t, fft_ubar, fft_vbar)

    fft_ubar .= champubar
    fft_vbar .= champvbar

    fft!(champu, 1)
    fft!(champv, 1)

    champu[1, :] .= 0.0
    champv[1, :] .= 0.0

    champu ./= (1im * m.ktau)
    champv ./= (1im * m.ktau)

    ifft!(champu, 1)
    ifft!(champv, 1)

    champu .*= m.epsilon
    champv .*= m.epsilon

    champu .+= transpose(champubar)
    champv .+= transpose(champvbar)

    champu .= ffu .- dtauh1u .- champu
    champv .= ffv .- dtauh1v .- champv

    fft!(champu, 1)
    fft!(champv, 1)

    fft_ug .= champu[1, :] ./ m.ntau
    fft_vg .= champv[1, :] ./ m.ntau

    champu[1, :] .= 0.0
    champv[1, :] .= 0.0

    champu ./= (1im .* m.ktau)
    champv ./= (1im .* m.ktau)

    ifft!(champu, 1)
    ifft!(champv, 1)

end

ext: jl


text:
DiffEqBase.has_reinit(i::AbstractNoiseProcess) = true
function DiffEqBase.reinit!(W::Union{NoiseProcess,NoiseApproximation},dt;
                            t0 = W.t[1],
                            erase_sol = true,
                            setup_next = false)

  if erase_sol
    resize!(W.t,1)
    resize!(W.W,1)
    if W.Z != nothing
      resize!(W.Z,1)
    end
  end

  W.curt = t0
  W.dt = dt
  if typeof(W) <: NoiseApproximation
    reinit!(W.source1)
    if W.source2 != nothing
      reinit!(W.source2)
    end
  end

  if isinplace(W)
    W.curW .= first(W.W)
    if W.Z != nothing
      W.curZ .= first(W.Z)
    end
  else
    W.curW = first(W.W)
    if W.Z != nothing
      W.curZ = first(W.Z)
    end
  end

  if typeof(W) <: NoiseProcess
    while !isempty(W.S₁)
      pop!(W.S₁) # Get a reset for this stack?
    end
    ResettableStacks.reset!(W.S₂)
  end
  setup_next && setup_next_step!(W)
  return nothing
end

function DiffEqBase.reinit!(W::AbstractNoiseProcess,dt;
                            t0 = W.t[1],
                            erase_sol = true,
                            setup_next = false)
  W.curt = t0
  W.dt = dt
  if typeof(W) <: NoiseGrid
    if isinplace(W)
      W.curW .= W.W[1]
      if W.Z != nothing
        W.curZ .= W.Z[1]
      end
    else
      W.curW = W.W[1]
      if W.Z != nothing
        W.curZ = W.Z[1]
      end
    end
    W.step_setup = true
  end
  setup_next && setup_next_step!(W)
  return nothing
end

function Base.reverse(W::AbstractNoiseProcess)
  if typeof(W) <: NoiseGrid
    backwardnoise = NoiseGrid(reverse(W.t),reverse(W.W))
  else
    W.save_everystep = false
    backwardnoise = NoiseWrapper(W, reverse=true)
  end
  return backwardnoise
end


text:
using PlotlyLight
using Test

@testset "history" begin
    PlotlyLight.set_history!(5)
    @test PlotlyLight.n_history[] == 5
    for _ in 1:6
        Plot(Config(; x = 1:10))
    end
    @test length(PlotlyLight.history()) == 5
end


text:
using TMLEEpistasis
using ArgParse

function parse_commandline()
    s = ArgParseSettings(description="Utilities to prepare genotypes for later PCA extraction.")

    @add_arg_table s begin
        "--input", "-i"
            help = "If merge is called, a common prefix to the bed files should be provided."*
                "If filter is called, a path to the chromosome .bed file should be provided"
            arg_type = String
            required = true
        "--output", "-o"
            help = "Path to the output file, if you use merge after filter you may want to use a common prefix"
            arg_type = String
            required = true
        "--qcfile"
            help = "Path to the UKBiobank ukb_snp_qc.txt"
            arg_type = String
        "--maf-threshold", "-t"
            help = "SNPs with MAF lower than this value will be filtered out"
            arg_type = Float64
            default = 0.01
            required = false
        "--ld-blocks", "-l"
            help = "Path to the LD block definition around the SNPs of interest. "*
                   "Issued from an LD analysis."
            arg_type = String
            required = false
        "filter"
            help = "Filter the SNPS from a PLINK .bed file based on quality metrics"
            action = :command
        "merge"
            help = "Merge multiple PLINK .bed files together"
            action = :command
        "adapt"
            help = "Adapts the csv output of flashpca to what is expected by the TMLE step"
            action = :command

    end

    return parse_args(s)
end

parsed_args = parse_commandline()

if haskey(parsed_args, "merge") == true
    merge_beds(parsed_args)
elseif haskey(parsed_args, "filter") == true
    filter_chromosome(parsed_args)
elseif haskey(parsed_args, "adapt")
    adapt_flashpca(parsed_args)
end


text:
const metre = 1e2
const meter = metre
const cm = metre / 1e2
const mm = metre / 1e3
const um = metre / 1e6
const nm = metre / 1e9

const second = 1e3
const ms = second / 1e3
const Hz = 1 / second

const voltage = 1e3
const mV = voltage / 1e3

const ampere = 1e6
const mA = ampere / 1e3
const uA = ampere / 1e6
const nA = ampere / 1e9

const farad=1e6
const uF = farad/1e6
const ufarad = uF

const siemens = 1e3
const mS = siemens / 1e3
const msiemens = mS
const nS = siemens / 1e9
const nsiemens = nS


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

"""
    OrdinaryKriging(γ)
    OrdinaryKriging(X, z, γ)

Ordinary Kriging with variogram model `γ`.

Optionally, pass the coordinates `X` and values `z`
to the [`fit`](@ref) function.
"""
struct OrdinaryKriging{G<:Variogram} <: KrigingEstimator
  γ::G
end

OrdinaryKriging(γ) = OrdinaryKriging{typeof(γ)}(γ)

OrdinaryKriging(X, z, γ) = GeoStatsBase.fit(OrdinaryKriging(γ), X, z)

nconstraints(estimator::OrdinaryKriging) = 1

function set_constraints_lhs!(estimator::OrdinaryKriging, LHS::AbstractMatrix, X::AbstractMatrix)
  T = eltype(LHS)
  LHS[end,:]   .= one(T)
  LHS[:,end]   .= one(T)
  LHS[end,end]  = zero(T)

  nothing
end

factorize(estimator::OrdinaryKriging, LHS::AbstractMatrix) = lu(LHS, check=false)

function set_constraints_rhs!(estimator::FittedKriging{E,S},
                              xₒ::AbstractVector) where {E<:OrdinaryKriging,S<:KrigingState}
  RHS = estimator.state.RHS

  RHS[end] = one(eltype(RHS))

  nothing
end


text:
# IMinuit initialization

const iminuit = PyNULL()
const mMinuit = PyNULL()
const iminuit_version = "iminuit=1.5.4"

# initialization -- anything that depends on Python has to go here,
# so that it occurs at runtime (while the rest can be precompiled).
# if PyCall is configured to use the Julia-specific Python, then iminuit
# can be automatically installed by pyimport_conda.
function __init__()
    copy!(iminuit, pyimport_conda("iminuit", iminuit_version, "conda-forge"))
    _version = iminuit.__version__
    if (_version < "1.5.0" || _version >= "1.6")
        println("The current iminuit version is " * _version * ". It will be changed to version 1.5.4.")
        run(`conda install iminuit=1.5.4 -c conda-forge`) # this updates the sys. conda, not the julia one.
    #     copy!(iminuit, pyimport_conda("iminuit", iminuit_version, "conda-forge"))
    # #     _version = iminuit.:version.__version__
    end
    # println("iminuit version " * _version * " has been imported as `iminuit`." )
    copy!(mMinuit, pyimport_conda("iminuit", iminuit_version, "conda-forge").:Minuit)
    # The following converts the ArgsView type to Vector{Float64}, but takes too much time
    # it takes more than 50 μs
    # pytype_mapping(minuit._libiminuit.ArgsView, Vector{Float64})
end


text:
using CountTimeSeries
using Test

@testset "CountTimeSeries.jl" begin
    using Random
    X = reshape(collect(1:100) .+ 0.0, (1, 100))
    X2 = vcat((1:100)', ((1:100).^2)') .+ 0.0
    models = Array{Any, 1}(undef, 30)

    models[1] = Model(pastObs = 1, pastMean = 1, X = X, external = [true], zi = true)
    models[2] = Model(pastObs = 1, X = X, external = [true], zi = true)
    models[3] = Model(model = "INARMA", pastObs = 1, pastMean = 1, X = X, external = [true], zi = true)
    models[4] = Model(model = "INARMA", pastObs = 1, X = X, external = [true], zi = true)
    models[5] = Model(model = "INARMA", pastMean = 1, X = X, external = [true], zi = true)
    models[6] = Model(pastObs = 1, pastMean = 1, distr = "NegativeBinomial", X = X, external = [true], zi = true)
    models[7] = Model(pastObs = 1, distr = "NegativeBinomial", X = X, external = [true], zi = true)
    models[8] = Model(model = "INARMA", pastObs = 1, pastMean = 1, distr = ["NegativeBinomial", "Poisson"], X = X, external = [true], zi = true)
    models[9] = Model(model = "INARMA", pastObs = 1, distr = "NegativeBinomial", X = X, external = [true], zi = true)
    models[10] = Model(model = "INARMA", pastMean = 1, distr = "NegativeBinomial", X = X, external = [true], zi = true)

    models[11] = Model(pastObs = 1, pastMean = 1, X = X2, external = [true, false], zi = true)
    models[12] = Model(pastObs = 1, X = X2, external = [true, false], zi = true)
    models[13] = Model(model = "INARMA", pastObs = 1, pastMean = 1, X = X2, external = [true, false], zi = true)
    models[14] = Model(model = "INARMA", pastObs = 1, X = X2, external = [true, false], zi = true)
    models[15] = Model(model = "INARMA", pastMean = 1, X = X2, external = [true, false], zi = true)
    models[16] = Model(pastObs = 1, pastMean = 1, distr = "NegativeBinomial", X = X2, external = [true, false], zi = true)
    models[17] = Model(pastObs = 1, distr = "NegativeBinomial", X = X2, external = [true, false], zi = true)
    models[18] = Model(model = "INARMA", pastObs = 1, pastMean = 1, distr = ["NegativeBinomial", "Poisson"], X = X2, external = [true, false], zi = true)
    models[19] = Model(model = "INARMA", pastObs = 1, distr = "NegativeBinomial", X = X2, external = [true, false], zi = true)
    models[20] = Model(model = "INARMA", pastMean = 1, distr = "NegativeBinomial", X = X2, external = [true, false], zi = true)

    models[21] = Model(pastObs = 1, pastMean = 1, X = X, external = [true], zi = true, link = "Log")
    models[22] = Model(pastObs = 1, X = X, external = [true], zi = true, link = "Log")
    models[23] = Model(model = "INARMA", pastObs = 1, pastMean = 1, X = X, external = [true], zi = true, link = "Log")
    models[24] = Model(model = "INARMA", pastObs = 1, X = X, external = [true], zi = true, link = "Log")
    models[25] = Model(model = "INARMA", pastMean = 1, X = X, external = [true], zi = true, link = "Log")
    models[26] = Model(pastObs = 1, pastMean = 1, distr = "NegativeBinomial", X = X, external = [true], zi = true, link = "Log")
    models[27] = Model(pastObs = 1, distr = "NegativeBinomial", X = X, external = [true], zi = true, link = "Log")
    models[28] = Model(model = "INARMA", pastObs = 1, pastMean = 1, distr = ["NegativeBinomial", "Poisson"], X = X, external = [true], zi = true, link = "Log")
    models[29] = Model(model = "INARMA", pastObs = 1, distr = "NegativeBinomial", X = X, external = [true], zi = true, link = "Log")
    models[30] = Model(model = "INARMA", pastMean = 1, distr = "NegativeBinomial", X = X, external = [true], zi = true, link = "Log")

    pars = Array{parameter, 1}(undef, 30)
    pars[1] = θ2par([10, 0.5, 0.2, 0.1, 0.1], models[1])
    pars[2] = θ2par([10, 0.5, 0.1, 0.1], models[2])
    pars[3] = θ2par([10, 0.5, 0.2, 0.1, 0.1], models[3])
    pars[4] = θ2par([10, 0.5, 0.1, 0.1], models[4])
    pars[5] = θ2par([10, 0.2, 0.1, 0.1], models[5])
    pars[6] = θ2par([10, 0.5, 0.2, 0.1, 3, 0.1], models[6])
    pars[7] = θ2par([10, 0.5, 0.1, 3, 0.1], models[7])
    pars[8] = θ2par([10, 0.5, 0.2, 0.1, 3, 0.1], models[8])
    pars[9] = θ2par([10, 0.2, 0.1, 3, 3, 0.1], models[9])
    pars[10] = θ2par([10, 0.5, 0.1, 3, 3, 0.1], models[10])

    pars[11] = θ2par([10, 0.5, 0.2, 0.1, 0.001, 0.1], models[11])
    pars[12] = θ2par([10, 0.5, 0.1, 0.001, 0.1], models[12])
    pars[13] = θ2par([10, 0.5, 0.2, 0.1, 0.001, 0.1], models[13])
    pars[14] = θ2par([10, 0.5, 0.1, 0.001, 0.1], models[14])
    pars[15] = θ2par([10, 0.2, 0.1, 0.001, 0.1], models[15])
    pars[16] = θ2par([10, 0.5, 0.2, 0.1, 0.001, 3, 0.1], models[16])
    pars[17] = θ2par([10, 0.5, 0.1, 0.001, 3, 0.1], models[17])
    pars[18] = θ2par([10, 0.5, 0.2, 0.1, 0.001, 3, 0.1], models[18])
    pars[19] = θ2par([10, 0.2, 0.1, 0.001, 3, 3, 0.1], models[19])
    pars[20] = θ2par([10, 0.5, 0.1, 0.001, 3, 3, 0.1], models[20])

    pars[21] = θ2par([1, 0.5, 0.2, -0.1, 0.1], models[21])
    pars[22] = θ2par([1, 0.5, -0.1, 0.1], models[22])
    pars[23] = θ2par([1, 0.5, 0.2, -0.1, 0.1], models[23])
    pars[24] = θ2par([1, 0.5, -0.1, 0.1], models[24])
    pars[25] = θ2par([1, 0.2, -0.1, 0.1], models[25])
    pars[26] = θ2par([1, 0.5, 0.2, -0.1, 3, 0.1], models[26])
    pars[27] = θ2par([1, 0.5, -0.1, 3, 0.1], models[27])
    pars[28] = θ2par([1, 0.5, 0.2, -0.1, 3, 0.1], models[28])
    pars[29] = θ2par([1, 0.2, -0.1, 3, 3, 0.1], models[29])
    pars[30] = θ2par([1, 0.5, -0.1, 3, 3, 0.1], models[30])

    Random.seed!(1)
    ys = (x -> x[1]).(simulate.(100, models, pars))
    LLs = (x -> x[1]).(ll.(ys, models, pars))

    settings = MLESettings.(ys, models, pars, ci = false)
    res = fit.(ys[[1, 2, 4, 6, 7, 9]], models[[1, 2, 4, 6, 7, 9]], settings[[1, 2, 4, 6, 7, 9]], printResults = true)

    res61 = fit(ys[6], models[1], printResults = false)
    res61qpois = QPois(res61)
    models[1] = Model(pastObs = 1, pastMean = 1, X = X, external = [true], zi = true)

    AIC.(res)
    BIC.(res)
    HQIC.(res)
    Xnew = reshape(collect(101:110) .+ 0.0, (1, 10))

    for i = 2:6
        predict(res[i], 10, 100, Xnew)
    end
    predict(res[2], 10, Xnew)
    pit(res[2])

    convert(INGARCHModel, models[2])
    convert(INGARCHModel, Model())
    convert(INARCHModel, Model())
    convert(INARMAModel, models[4])

    show.(models)
    show.(res)
    show.(pars)

    model = Model(model = "INARMA", pastMean = 1:2)
    y = simulate(100, model, [10, 0.5, 0.2])[1]
    ll(y, model, [10, 0.5, 0.2])

    model = Model(model = "INARMA", pastMean = 1:2, pastObs = 1, zi = true)
    y = simulate(100, model, [10, 0.1, 0.5, 0.2, 0.1])[1]
    ll(y, model, [10, 0.1, 0.5, 0.2, 0.1])

    model = Model(zi = true)
    y = simulate(100, model, [10.0, 0.1])[1]
    ll(y, model, [10, 0.1])

    model = Model(zi = true, X = X2, external = [true, true])
    y = simulate(100, model, [10, 0.1, 0.01, 0.01])[1]
    ll(y, model, [10, 0.1, 0.01, 0.01])
end


text:
function hm_create_integrationvolumes!(
    result::HMIData{T, I, HyperRectVolume{T}},
    settings::HMISettings) where {T<:AbstractFloat, I<:Integer}

    if isempty(result.volumelist1) || isempty(result.volumelist2)
        nvols = (isempty(result.volumelist1) ? length(result.dataset1.startingIDs) : 0) +
                (isempty(result.volumelist2) ? length(result.dataset2.startingIDs) : 0)

        @info "Create $nvols Hyperrectangles using $(_global_mt_setting ? nthreads() : 1) thread(s)"
        progressbar = Progress(nvols)

        isempty(result.volumelist1) && hm_create_integrationvolumes_dataset!(
            result.dataset1, result.volumelist1, result.cubelist1, result.iterations1, result.whiteningresult.targetprobfactor,
                progressbar, settings)
        isempty(result.volumelist2) && hm_create_integrationvolumes_dataset!(
            result.dataset2, result.volumelist2, result.cubelist2, result.iterations2, result.whiteningresult.targetprobfactor,
                progressbar, settings)

        finish!(progressbar)
    end

    nvols = (result.dataset1.isnew ? length(result.volumelist1) : 0) +
            (result.dataset2.isnew ? length(result.volumelist2) : 0)

    if nvols > 0
        @info "Updating $nvols Hyperrectangles using $(_global_mt_setting ? nthreads() : 1) thread(s)"
        progressbar = Progress(nvols)

        result.dataset2.isnew && hm_update_integrationvolumes_dataset!(result.dataset2, result.volumelist1, progressbar)
        result.dataset1.isnew && hm_update_integrationvolumes_dataset!(result.dataset1, result.volumelist2, progressbar)

        finish!(progressbar)
    end
end

function hm_create_integrationvolumes_dataset!(
    dataset::DataSet{T, I},
    volumes::Array{IntegrationVolume{T, I, HyperRectVolume{T}}, 1},
    cubes::Array{HyperRectVolume{T}, 1},
    iterations_per_volume::Array{I, 1},
    targetprobfactor::T,
    progressbar::Progress,
    settings::HMISettings) where {T<:AbstractFloat, I<:Integer}

    maxPoints::I = 0
    totalpoints::I = 0

    thread_volumes = Vector{IntegrationVolume{T, I, HyperRectVolume{T}}}(undef, length(dataset.startingIDs))
    thread_cubes = Vector{HyperRectVolume{T}}(undef, length(dataset.startingIDs))
    thread_iterations = Vector{I}(undef, length(dataset.startingIDs))

    atomic_centerID = Atomic{I}(1)

    @mt BAT.hyperrectangle_creationproccess!(dataset, targetprobfactor, settings,
        thread_volumes, thread_cubes, thread_iterations, atomic_centerID, progressbar)

    for i in eachindex(thread_volumes)
        if isassigned(thread_volumes, i) == false
        elseif thread_volumes[i].pointcloud.probfactor == 1.0 || thread_volumes[i].pointcloud.points < dataset.P * 4
        else
            push!(volumes, thread_volumes[i])
            push!(cubes, thread_cubes[i])
            push!(iterations_per_volume, thread_iterations[i])
            maxPoints = max(maxPoints, thread_volumes[i].pointcloud.points)
            totalpoints += thread_volumes[i].pointcloud.points
        end
    end

    dataset.isnew = true
end

function hm_update_integrationvolumes_dataset!(
    dataset::DataSet{T, I},
    volumes::Array{IntegrationVolume{T, I, HyperRectVolume{T}}, 1},
    progressbar::Progress) where {T<:AbstractFloat, I<:Integer}

    maxPoints = zero(T)

    wp_volumes = workpart(volumes, ParallelProcessingTools.workers(), ParallelProcessingTools.myid())
    @mt for i in workpart(eachindex(wp_volumes), mt_nthreads(), threadid())
        update!(volumes[i], dataset)

        maxPoints = max(maxPoints, volumes[i].pointcloud.points)

        @critical next!(progressbar)
    end

    #remove rectangles with less than 1% points of the largest rectangle (in terms of points)
    j = length(volumes)
    for i = 1:length(volumes)
        if volumes[j].pointcloud.points < maxPoints * 0.01 || volumes[j].pointcloud.points < dataset.P * 4
            deleteat!(volumes, j)
        end
        j -= 1
    end
    dataset.isnew = false
end

"""
    hyperrectangle_creationproccess!(...)

*AHMI-internal, not part of stable public API.*

This function assigns each thread its own hyper-rectangle to build, if in multithreading-mode.
"""
function hyperrectangle_creationproccess!(
    dataset::DataSet{T, I},
    targetprobfactor::T,
    settings::HMISettings,
    integrationvolumes::Vector{IntegrationVolume{T, I, HyperRectVolume{T}}},
    cubevolumes::Vector{HyperRectVolume{T}},
    total_iterations::Vector{I},
    atomic_centerID::Atomic{I},
    progressbar::Progress) where {T<:AbstractFloat, I<:Integer}

    while true
        #get new starting sample
        idc = atomic_add!(atomic_centerID, 1)
        if idc > length(dataset.startingIDs)
            break
        end
        id = dataset.startingIDs[idc]

        #update progress bar
        @critical next!(progressbar)

        #create hyper-rectangle
        integrationvolumes[idc], cubevolumes[idc], total_iterations[idc] = create_hyperrectangle(id, dataset, targetprobfactor, settings)

        @debug "Hyperrectangle created. Points:\t$(integrationvolumes[idc].pointcloud.points)\tVolume:\t$(integrationvolumes[idc].volume)\tProb. Factor:\t$(integrationvolumes[idc].pointcloud.probfactor)"
    end
end


text:
using Genie, Genie.Router
using Genie.Renderer, Genie.Renderer.Html, Genie.Renderer.Json
using JSON
using SwagUI

swagger_document = JSON.parsefile("./swagger.json")

urls = Array{Dict{String, Any}, 1}()
url1 = Dict{String, Any}()
url1["url"] = "https://petstore.swagger.io/v2/swagger.json"
url1["name"] = "Spec1"
url2 = Dict{String, Any}()
url2["url"] = "https://petstore.swagger.io/v2/swagger.json"
url2["name"] = "Spec2"
push!(urls, url1)
push!(urls, url2)

options = Options()
# options.custom_css = ".swagger-ui .topbar { display: none }"
options.show_explorer = true

# options.swagger_options["validatorUrl"] = nothing
options.swagger_options["url"] = "https://petstore.swagger.io/v2/swagger.json"
# options.swagger_options["urls"] = urls

route("/docs") do
    render_swagger(nothing, options=options)
    # render_swagger(swagger_document, options=options)
end

up(8001, async = false)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    TransformedDensity <: AbstractDensity

Abstract type for transformed densities.

In addition to the [`AbstractDensity`](@ref) interface, subtypes must support
the functions `Base.parent` and [`trafoof`](@ref):

```julia
parent(d::SomeTransformedDensity)::AbstractDensity
trafoof(d::SomeTransformedDensity)::Function
```
"""
struct AbstractTransformedDensity end
export AbstractTransformedDensity

"""
    trafoof(d::AbstractTransformedDensity)::AbstractDensity

Get the transform from `parent(d)` to `d`, so that

```julia
trafoof(d)(parent(d)) == d
```
"""
function trafoof end
export trafoof

abstract type TDVolCorr end
struct TDNoCorr <: TDVolCorr end
struct TDLADJCorr <: TDVolCorr end

"""
    TransformedDensity

*BAT-internal, not part of stable public API.*
"""
struct TransformedDensity{D<:AbstractDensity,FT<:Function,VC<:TDVolCorr,VS<:AbstractValueShape} <: AbstractDensity
    orig::D
    trafo::FT  # ToDo: store inverse(trafo) instead?
    volcorr::VC
    _varshape::VS
end

function TransformedDensity(orig::AbstractDensity, trafo::Function, volcorr::TDVolCorr)
    vs = trafo(varshape(orig))
    TransformedDensity(orig, trafo, volcorr, vs)
end

@inline function (trafo::DistributionTransform)(density::AbstractDensity; volcorr::Val{vc} = Val(true)) where vc
    if vc
        TransformedDensity(density, trafo, TDLADJCorr())
    else
        TransformedDensity(density, trafo, TDNoCorr())
    end
end

Base.parent(density::TransformedDensity) = density.orig
trafoof(density::TransformedDensity) = density.trafo

ValueShapes.varshape(density::TransformedDensity) = density._varshape

# ToDo: Should not be neccessary, improve default implementation of
# ValueShapes.totalndof(density::AbstractDensity):
ValueShapes.totalndof(density::TransformedDensity) = totalndof(varshape(density))

var_bounds(density::TransformedDensity{<:Any,<:DistributionTransform}) = dist_param_bounds(density.trafo.target_dist)

function DensityInterface.logdensityof(density::TransformedDensity{D,FT,TDNoCorr}, v::Any) where {D,FT}
    v_orig = inverse(density.trafo)(v)
    logdensityof(parent(density), v_orig)
end

function checked_logdensityof(density::TransformedDensity{D,FT,TDNoCorr}, v::Any) where {D,FT}
    v_orig = inverse(density.trafo)(v)
    checked_logdensityof(parent(density), v_orig)
end

function _v_orig_and_ladj(density::TransformedDensity, v::Any)
    with_logabsdet_jacobian(inverse(density.trafo), v)
end

# TODO: Would profit from custom pullback:
function _combine_logd_with_ladj(logd_orig::Real, ladj::Real)
    logd_result = logd_orig + ladj
    R = typeof(logd_result)

    if isnan(logd_result) && logd_orig == -Inf && ladj == +Inf
        # Zero density wins against infinite volume:
        R(-Inf)
    elseif isfinite(logd_orig) && (ladj == -Inf)
        # Maybe  also for (logd_orig == -Inf) && isfinite(ladj) ?
        # Return constant -Inf to prevent problems with ForwardDiff:
        #R(-Inf)
        near_neg_inf(R) # Avoids AdvancedHMC warnings
    else
        logd_result
    end
end

function DensityInterface.logdensityof(density::TransformedDensity{D,FT,TDLADJCorr}, v::Any) where {D,FT,}
    v_orig, ladj = _v_orig_and_ladj(density, v)
    logd_orig = logdensityof(parent(density), v_orig)
    _combine_logd_with_ladj(logd_orig, ladj)
end

function checked_logdensityof(density::TransformedDensity{D,FT,TDLADJCorr}, v::Any) where {D,FT,}
    v_orig, ladj = _v_orig_and_ladj(density, v)
    logd_orig = checked_logdensityof(parent(density), v_orig)
    _combine_logd_with_ladj(logd_orig, ladj)
end


text:
function determine_impedance_parameters(input_data, spatial_data, spatial_data_matrices, optimal_path, equipment_data)
    r = 0
    x = 0
    bc = 0
    r_pu = 0
    x_pu = 0
    bc_pu = 0
    km = 0
    km_ohl = 0
    km_ugc = 0
    for idx = 1 : size(optimal_path, 1) - 1
        ohl_cable = ohl_or_cable(spatial_data, spatial_data_matrices, optimal_path[idx], optimal_path[idx + 1])
        r_, x_, bc_, r_pu_, x_pu_, bc_pu_ = get_impedance_data(equipment_data, ohl_cable, input_data)
        km_ = get_segment_distance(spatial_data, spatial_data_matrices, optimal_path[idx], optimal_path[idx + 1])

        km = km + km_
        if ohl_cable == "ohl"
            km_ohl = km_ohl + km_
        else
            km_ugc = km_ugc + km_
        end
        r = r + (r_ * km_)
        r_pu = r_pu + (r_pu_ * km_)
        x = x + (x_ * km_)
        x_pu = x_pu +(x_pu_ * km_)

        bc = bc + (bc_ * km_)
        bc_pu = bc_pu + (bc_pu_ * km_)
    end

    route_impedance = Dict{String, Any}()
    route_length = Dict{String, Any}()

    route_length["total_length"] = km
    route_length["ohl_length"] = km_ohl
    route_length["ugc_length"] = km_ugc

    route_impedance["r"] = r
    route_impedance["x"] = x
    route_impedance["bc"] = bc

    route_impedance["r_pu"] = r_pu
    route_impedance["x_pu"] = x_pu
    route_impedance["bc_pu"] = bc_pu

    return route_impedance, route_length
 end

 function get_impedance_data(equipment_data::Dict, ohl_cable::String, input_data::Dict)
    technology = input_data["technology"]
    cross_section = equipment_data[ohl_cable]["onshore"]["cross_section"]
    circuits = equipment_data[ohl_cable]["onshore"]["circuits"]
    if ohl_cable == "ohl"
        bundles = input_data[ohl_cable]["onshore"]["bundles"]
    else
        bundles = 1
    end
    r = input_data["impedances"][technology][ohl_cable]["$cross_section"]["r"] / (bundles * circuits)
    x = input_data["impedances"][technology][ohl_cable]["$cross_section"]["x"] / (bundles * circuits)
    bc = input_data["impedances"][technology][ohl_cable]["$cross_section"]["c"] * circuits * 10e-9 * 2 * pi * 50

    v_str = join([technology, "_", ohl_cable])
    voltage = input_data["voltages"][v_str]
    Zbase = voltage^2 / 100

    r_pu = r / Zbase
    x_pu = x / Zbase
    bc_pu = bc * Zbase

    return r, x, bc, r_pu, x_pu, bc_pu
 end

 function get_segment_distance(spatial_data::Dict, spatial_data_matrices::Dict, node_id1::Int, node_id2::Int)
    index_s1 = findall(node_id1 .== spatial_data["segments"][:, 2])
    index_s2 = findall(node_id2 .== spatial_data["segments"][:, 3])
    index_s = intersect(index_s1, index_s2)
    if isempty(index_s)
        index_s1 = findall(node_id2 .== spatial_data["segments"][:, 2])
        index_s2 = findall(node_id1 .== spatial_data["segments"][:, 3])
        index_s = intersect(index_s1, index_s2)
    end

    if !isempty(findall(index_s .== spatial_data_matrices["ohl"]["segments"][:, 1]))
        km = spatial_data_matrices["ohl"]["segment_km"][index_s]
    else
        index_ugc = index_s[1] - size(spatial_data_matrices["ohl"]["segments"], 1)
        if !isempty(findall(index_s .== spatial_data_matrices["ugc"]["segments"][:, 1]))
            km = spatial_data_matrices["ugc"]["segment_km"][index_ugc]
        else
            km = 0
        end
    end

    return km
 end


text:
include("../roblib.jl")

function get_average_temperature(r)
    count = 0
    all_temperature = 0
    direction = Right

    while ((isborder(r, Right) == false) || (isborder(r, Up) == false))
        if ismarker(r)
            count += 1
            all_temperature += temperature(r)
        end
        if isborder(r, direction)
            direction = inverse(direction)
            move!(r, Up)
        else move!(r, direction)
        end
    end

    println(all_temperature/count)
end


text:
import StatsBase: sample, WeightVec

export SBM, sample

type SBM <: Graphon
    alpha::Array{Float64,1}
    P::Array{Float64,2}
#     function SBM(alpha::Array{Float64,1}, P::Array{Float64,2})
#         assert(size(P,1)==size(P,2))
#         assert(P==P')
#         assert(length(alpha)==size(P,1))
#         assert(max(P)<=1)
#         assert(min(P)>=0)
#         assert(sum(alpha)==1)
#         assert(alpha>=0)
#         return new(alpha, P)
#     end
end
function SBM(n::Int)
    alpha = rand(n)
    alpha /= sum(alpha)
    P = symmetrize!(rand(n,n))
    return SBM(alpha, P)
end
function SBM(alpha::Array{Float64,1})
    alpha /= sum(alpha)
    n = length(alpha)
    P = symmetrize!(rand(n,n))
    return SBM(alpha,P)
end
function SBM(P::Array{Float64,2})
    n = size(P,1)
    alpha = ones(n)/n
    return SBM(alpha,P)
end

## Sampling

function sample(sbm::SBM, nsamples::Int)
    n = length(sbm.alpha)
    indices = sample([1:n], WeightVec(sbm.alpha), nsamples)
    R = symmetrize!(rand(nsamples, nsamples))
    A = int(R .<= sbm.P[indices, indices])
    # remove self loops
    for i=1:n
        A[i,i] = 0
    end
    return A
end


text:
### Routines for constructing smooth functions that generate field data ###

import Base: *, -, +, show

export AbstractSpatialField, Gaussian, radius, center, strength,
       EmptySpatialField,
       SpatialGaussian,
       GeneratedField, datatype, grid,
       PulseField

abstract type AbstractSpatialField end
abstract type AbstractGeneratedField end

## Empty spatial field

"""
    EmptySpatialField()

Create a blank spatial field. This is primarily useful for initializing a
sum of spatial fields.

# Example
```jldoctest
julia> g = EmptySpatialField()
EmptySpatialField()

julia> g(2,3)
0.0
```
"""
struct EmptySpatialField <: AbstractSpatialField end
(g::EmptySpatialField)(x,y) = Float64(0)

## Gaussian

"""
    Gaussian(σ,x0,A)

Construct a 1-d Gaussian function centered at `x0` with standard deviation `σ`
and amplitude `A`. The resulting function can be evaluated at any real-valued
number.

# Example

```jldoctest
julia> g = Gaussian(0.2,0,1)
Gaussian(0.2, 0, 1, 2.8209479177387813)

julia> g(0.2)
1.0377687435514866
```
"""
struct Gaussian
  σ :: Real
  x0 :: Real
  A :: Real
  fact :: Float64
end

Gaussian(σ,x0,A) = Gaussian(σ,x0,A,A/sqrt(π*σ^2))

radius(g::Gaussian) = g.σ
center(g::Gaussian) = g.x0
strength(g::Gaussian) = g.A

@inline gaussian(r;tol=6.0) = abs(r) < tol ? exp(-r^2) : 0.0

(g::Gaussian)(x) = g.fact*gaussian((x-center(g))/radius(g))

## Spatial Gaussian field ##

struct SpatialGaussian <: AbstractSpatialField
  gx :: Gaussian
  gy :: Gaussian
  A :: Real
end

SpatialGaussian(σx::Real,σy::Real,x0::Real,y0::Real,A::Real) =
          SpatialGaussian(Gaussian(σx,x0,A),Gaussian(σy,y0,1),A)

SpatialGaussian(σ::Real,x0::Real,y0::Real,A::Real) = SpatialGaussian(σ,σ,x0,y0,A)

(g::SpatialGaussian)(x,y) = g.gx(x)*g.gy(y)

## Scaling spatial fields

struct ScaledField{N <: Real, P <: AbstractSpatialField} <: AbstractSpatialField
    s::N
    p::P
end
function show(io::IO, p::ScaledField)
    print(io, "$(p.s) × ($(p.p))")
end
s::Number * p::AbstractSpatialField = ScaledField(s, p)
-(p::AbstractSpatialField) = ScaledField(-1, p)
(p::ScaledField)(x,y) = p.s*p.p(x,y)

struct AddedFields{T <: Tuple} <: AbstractSpatialField
    ps::T
end
function show(io::IO, Σp::AddedFields)
    println(io, "AddedFields:")
    for p in Σp.ps
        println(io, "  $p")
    end
end

## Adding spatial fields together

"""
    p₁::AbstractSpatialField + p₂::AbstractSpatialField

Add the fields so that `(p₁ + p₂)(x,y) = p₁(x,y) + p₂(x,y)`.

"""
+(p::AbstractSpatialField, Σp::AddedFields) = AddedFields((p, Σp.ps...))
+(Σp::AddedFields, p::AbstractSpatialField) = AddedFields((Σp.ps..., p))
function +(Σp₁::AddedFields, Σp₂::AddedFields)
    AddedFields((Σp₁..., Σp₂...))
end

-(p₁::AbstractSpatialField, p₂::AbstractSpatialField) = p₁ + (-p₂)

+(p::AbstractSpatialField...) = AddedFields(p)

function (Σp::AddedFields)(x,y)
    f = 0.0
    for p in Σp.ps
        f += p(x,y)
    end
    f
end

## For generating an actual instance of a spatial field
"""
    GeneratedField(d::GridData,field::AbstractSpatialField...,grid::PhysicalGrid)

Create an instance of a spatial field function `field` on scalar grid data `d`,
based on a grid `grid`. After creating the instance `g = GeneratedField(d,field,grid)`,
then the resulting grid data can be accessed by typing `g()`. For vector grid data,
a separate `field` must be supplied for each component.
"""
struct GeneratedField{T <: GridData}
    fielddata :: T
    fieldfcns :: Vector{AbstractSpatialField}
    grid :: PhysicalGrid
end

function GeneratedField(d::ScalarGridData,field::AbstractSpatialField,g::PhysicalGrid)
    xg, yg = coordinates(d,g)
    tmp = _generatedfield(xg,yg,field,d)
    GeneratedField{typeof(tmp)}(tmp,AbstractSpatialField[field],g)
end

function GeneratedField(d::CollectedGridData,fields::Vector{AbstractSpatialField},g::PhysicalGrid)
  @assert length(fields) == _numberofcomponents(typeof(d))
  tmp = similar(d)
  cnt = 0
  for (i,fname) in enumerate(propertynames(d))
    if typeof(getfield(d,fname)) <: GridData
        cnt += 1
        x, y = coordinates(getfield(d,fname),g)
        getfield(tmp,fname) .= _generatedfield(x,y,fields[cnt],getfield(d,fname))
    end
  end
  GeneratedField{typeof(tmp)}(tmp,fields,g)
end

GeneratedField(d::VectorGridData,
      fieldu::AbstractSpatialField,fieldv::AbstractSpatialField,
      g::PhysicalGrid) = GeneratedField(d,AbstractSpatialField[fieldu,fieldv],g)

_generatedfield(xg,yg,field::AbstractSpatialField,d::ScalarGridData) =
          typeof(d)(field.(xg*ones(1,length(yg)),ones(length(xg))*yg'))

datatype(f::GeneratedField{T}) where {T} = T
grid(f::GeneratedField{T}) where {T} = f.grid

(f::GeneratedField{T})() where {T} = f.fielddata

## For generating a transient form of a spatial field

"""
    PulseField(g::GeneratedField,t0,σ)

Create a transient form of a generated spatial field, useful for
introducing a short-lived forcing field onto the grid. The supplied field `g`
is modulated by a Gaussian pulse centered at time `t0` with width `σ`. This pulse's
maximum is unity, so that the amplitude of the overall field is set by `g`. The
resulting object can be evaluated with a single argument (time) and returns a
`ScalarGridData` type object of the same type contained in `g`.
"""
struct PulseField
  gfield :: GeneratedField
  timemod :: Gaussian
end

PulseField(gfield::GeneratedField{T},t0::Real,σt::Real) where {T <: GridData} =
            PulseField(gfield,Gaussian(σt,t0,sqrt(π*σt^2)))

(f::PulseField)(t) = f.timemod(t)*f.gfield()

datatype(f::PulseField) = datatype(f.gfield)

## Allow all spatial fields to accept a vector argument

#for T in InteractiveUtils.subtypes(AbstractSpatialField)
#    @eval (f::$T)(x::Vector{R}) where {R <: Float64} = f(x...)
#end


text:
module Common

using ArgParse
using JSON

export parse_commandline
export SEGMENT_ID, MIP, VOXEL_SIZE, SEGMENTATION_LAYER
const MIP = UInt32(3)
const VOXEL_SIZE = (5,5,45)

function __init__()
    # setup AWS secrets
    if isfile("/secrets/aws-secret.json")
        d = JSON.parsefile("/secrets/aws-secret.json")
        for (k,v) in d
            ENV[k] = v
        end
    end
end

"""
customized argument parse for tuple
"""
function ArgParse.parse_item(::Type{NTuple{3,Int}}, x::AbstractString)
    return (map(Meta.parse, split(x,","))...,)
end
function ArgParse.parse_item(::Type{UInt32}, x::AbstractString)
    return UInt32(Meta.parse(x))
end

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--neuronid", "-i"
            help = "the segment id to skeletonize"
            arg_type = Int
            # default = SEGMENT_ID #77497
        "--swcdir", "-s"
            help = "the directory to store swc file"
            arg_type = String
            default = "/tmp"
        "--outputpath", "-o"
            help = "the path to store outputs"
            arg_type = String
            default = "/tmp"
        "--voxelsize", "-v"
            help = "voxel size of the raw image, mip level 0"
            arg_type = NTuple{3,Int}
            default = VOXEL_SIZE
        "--mip", "-m"
            help = "mip level of the dataset"
            arg_type = UInt32
            default = MIP
        "--meshname", "-e"
            help = "folder name of mesh files"
            arg_type = String
            default = "mesh_mip_4"
        "--idlistfile", "-f"
            help = "the id list in text from google spreasheet"
            arg_type = String
        "--sqsqueue", "-q"
            help = "AWS SQS queue name"
            arg_type = String
        "--segmentationlayer", "-l"
            help = "segmentation layer path in the cloud storage, only support Google Cloud Storage now"
            arg_type = String
    end
    return parse_args(s)
end

end # module


text:
# generic lattice plotting
function plot(
    		lattice :: L
            ;
            site_labels :: Bool = false,
    		site_radius :: Real = 25,
    		bond_thickness :: Real = 8,
    		visualize_periodic :: Bool = false,
    		colorcode_sites :: Union{Symbol,Dict} = Dict(),
    		colorcode_bonds :: Union{Symbol,Dict} = Dict(),
            kwargs...
        ) where {LS,D,LB,N,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N},U,L<:AbstractLattice{S,B,U}}

    ##########-------------------------------
    # STEP 1 #  set the colors correctly
    ##########-------------------------------

    # set the fallback colors (site)
    color_fallback_site = Int64[50,50,50]
    label_fallback_site = getDefaultLabel(LS)

    # set the fallback colors (bond)
    color_fallback_bond = Int64[10,10,10]
    label_fallback_bond = getDefaultLabel(LB)

    # TODO automatically set dictonaries
    if typeof(colorcode_bonds) == Symbol
        if colorcode_bonds == :kitaev || colorcode_bonds == :Kitaev
            colorcode_bonds = generateBondColorcodeKitaev(lattice)
        else
            println("colorcode :$(colorcode_bonds) could not be built, using fallback")
            colorcode_bonds = Dict()
        end
    end
    if typeof(colorcode_sites) == Symbol
        println("colorcode :$(colorcode_sites) could not be built, using fallback")
        colorcode_sites = Dict()
    end

	# repair color dictonaries with fallback entries
	#colorcode_bonds[label_fallback_bond] = get(colorcode_bonds, label_fallback_bond, color_fallback_bond)
	#colorcode_sites[label_fallback_site] = get(colorcode_sites, label_fallback_site, color_fallback_site)

    ##########-------------------------------
    # STEP 2 #  Plotting
    ##########-------------------------------

    ############---------------------
    # STEP 2.1 #  Plot all Bonds
    ############---------------------

    # get all bond labels
    bond_label_list = unique!(label.(bonds(lattice)))
    # iterate over all bond labels
    for l in bond_label_list
        # get all bonds to that label
        bond_list = B[b for b in bonds(lattice) if label(b)==l && from(b)>to(b) && (!isPeriodic(b) || visualize_periodic)]
        color     = get(colorcode_bonds, l, color_fallback_bond)
        # plot the bonds
        plotBonds(
            bond_list,
            sites(lattice),
            bond_thickness,
            color;
            kwargs...
        )
    end

    ############---------------------
    # STEP 2.2 #  Plot all Sites
    ############---------------------

    # get all site labels
    site_label_list = unique!(label.(sites(lattice)))

    # iterate over all site labels
    for l in site_label_list
        # get all sites to that label
        site_list = S[s for s in sites(lattice) if label(s)==l]
        color     = get(colorcode_sites, l, color_fallback_site)
        # plot the site
        plotSites(
            site_list,
            site_radius,
            color;
            site_labels = site_labels,
            kwargs...
        )
    end

end

# export the plot function
export plot


text:
#==
Example taken from
Messac et. al., Amir Ismail-Yahaya, and Christopher A Mattson.
The normalized normal constraint method for generating the pareto frontier.
Structural and multidisciplinary optimization, 25(2):86–98, 2003.

This example shows that the NBI method will generate local Pareto points that
are not global and dominated points that are not locally optimal.
We can make sure all NBI points are locally Pareto optimal by using an
inequality-constrained extension.
==#

using MultiJuMP, JuMP
using Ipopt

m = MultiModel(solver = IpoptSolver())
@variable(m, 0 <= x[i=1:2] <= 5)
@NLexpression(m, f1, x[1])
@NLexpression(m, f2, x[2])
@NLconstraint(m, 5exp(-x[1])+2exp(-0.5(x[1]-3)^2) <= x[2])

obj1 = SingleObjective(f1, sense = :Min)

# As the problem is nonconvex, we have to supply the
# initial value to get the global minimum of f2
obj2 = SingleObjective(f2, sense = :Min,
                       iv = Dict{Symbol, Any}(:x => [5., 0.]))

multim = getMultiData(m)
multim.objectives = [obj1, obj2]
multim.pointsperdim = 60

solve(m, method = :NBI, inequalityconstraint = true)
#solve(m, method = :NBI)

function pltfun(x1)
    setvalue(x[1], x1)
    setvalue(x[2:n], zeros(n-1))
    getvalue(f2)
end
pltfun(x1) = 5exp(-x1)+2exp(-0.5(x1-3)^2)

plt = plot(multim, label="NBI points")
plot!(plt, pltfun, 0,5, markercolor="orange",
      label="Image boundary")


text:
include("../src/HelloDB.jl")

module TestHelpers
  using HelloDB: setdb, dropdb, DB

  function set_test_db()
    setdb("test_db")
  end

  function drop_test_db()
    if DB != false
      dropdb(DB)
    end
  end

  function src_require(filename)
    require("../src/$(filename)")
  end
end


text:
"""
	grdcut(cmd0::String="", arg1=[], kwargs...)

Produce a new outgrid which is a subregion of ingrid. The subregion is specified with
``limits`` (the -R); the specified range must not exceed the range of ingrid (but see ``extend``).

Full option list at [`grdcut`]($(GMTdoc)grdcut.html)

Parameters
----------

- **F** | **clip** | **cutline** :: [Type => Str | GMTdaset | Mx2 array | NamedTuple]	`Arg = array|fname[+c] | (polygon=Array|Str, crop2cutline=Bool, invert=Bool)`

    Specify a closed polygon (either a file or a dataset). All grid nodes outside the
    polygon will be set to NaN (>= GMT6.2).
    ($(GMTdoc)grdcut.html#f)
- **G** | **save** | **outgrid** | **outfile** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdcut(....) form.
    ($(GMTdoc)grdcut.html#g)
- **img** | **usegdal** | **gdal** :: [Type => Any]

    Force the cut operation to be done by GDAL. Works for images where GMT fails or even crash.
- $(GMT.opt_J)
- **N** | **extend** :: [Type => Str or []]

    Allow grid to be extended if new region exceeds existing boundaries. Append nodata value
    to initialize nodes outside current region [Default is NaN].
    ($(GMTdoc)grdcut.html#n)
- $(GMT.opt_R)
- **S** | **circ_subregion** :: [Type => Str]    ``Arg = [n]lon/lat/radius[unit]``

    Specify an origin and radius; append a distance unit and we determine the corresponding
    rectangular region so that all grid nodes on or inside the circle are contained in the subset.
    ($(GMTdoc)grdcut.html#s)
- $(GMT.opt_V)
- **Z** | **z_subregion** :: [Type => Str]       ``Arg = [n|N |r][min/max]``

    Determine a new rectangular region so that all nodes outside this region are also outside
    the given z-range.
    ($(GMTdoc)grdcut.html#z)
- $(GMT.opt_f)
"""
function grdcut(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && contains(cmd0, " ") && return monolitic("grdcut", cmd0, arg1)

	arg2 = nothing
	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode
	cmd, opt_R = parse_R(d, "")
	(opt_R == "") && error("Must provide the cutting limits (GMT option R)")
    cmd, = parse_common_opts(d, cmd, [:V_params :f])
    opt_J, = parse_J(d, "")
    (!startswith(opt_J, " -JX")) && (cmd *= opt_J)
	cmd = parse_these_opts(cmd, d, [[:D], [:N :extend], [:S :circ_subregion], [:Z :z_subregion]])
	opt_G = parse_G(d, "")[1]
	outname = (opt_G != "") ? opt_G[4:end] : ""
	cmd *= opt_G
	cmd, args, n, = add_opt(d, cmd, "F", [:F :clip :cutline], :polygon, Array{Any,1}([arg1, arg2]),
	                        (crop2cutline="_+c", invert="_+i"))
	if (n > 0)  arg1, arg2 = args[:]  end
	(show_kwargs[1]) && return print_kwarg_opts([:img :usegdal :gdal], "Any")		# Just print the options

	if (cmd0 != "" && (guess_T_from_ext(cmd0) == " -Ti" || (find_in_dict(d, [:usegdal :gdal])[1]) !== nothing))
		(dbg_print_cmd(d, cmd) !== nothing) && return "grdcut $cmd0 " * cmd		# Vd=2 cause this return
		t = split(scan_opt(cmd, "-R"), '/')
		opts = ["-projwin", t[1], t[4], t[2], t[3]]		# -projwin <ulx> <uly> <lrx> <lry>
		cut_with_gdal(cmd0, opts, outname)
	else
		common_grd(d, cmd0, cmd, "grdcut ", arg1, arg2)	# Finish build cmd and run it
	end
end

function cut_with_gdal(fname::String, opts::Vector{<:AbstractString}, outname::String="")
	if (outname == "")
		gdaltranslate(fname, opts)	# Layout is "TRB" so all matrices are contrary to Julia order
	else
		gdaltranslate(fname, opts; dest=outname)
		return nothing				# Since it wrote a file so nothing to return
	end
end

# ---------------------------------------------------------------------------------------------------
grdcut(arg1, cmd0::String=""; kw...) = grdcut(cmd0, arg1; kw...)

# ---------------------------------------------------------------------------------------------------
"""
crop(arg::GItype; kw...)

Crop a subregion of a grid (GMTgrid) or a image (GMTimage). The subregion is specified with the
``limits`` or ``region`` keyword; the specified range must not exceed the range of the input.
This function differs from ``grdcut`` in the sense that it doesn't call the GMT lib and works only on
in-memory array (i.e., no disk files).

### Returns
A grid or an image, depending on the input type, plus two 1x2 matrices with the indices of the cropped zone.

## Example
	G = GMT.peaks();
	crop(G, region=(-2,2,-2,2))
"""
function crop(arg::GItype; kw...)
	d = KW(kw)
	_, opt_R = parse_R(d, "")
	(opt_R == "") && error("Must provide the cropping limits")
	lims = opt_R2num(opt_R)
	# Must test that requested cropping limits fit inside array BB
	lims[1], lims[2] = max(lims[1], arg.range[1]), min(lims[2], arg.range[2])	# Avoid overflows in Region
	lims[3], lims[4] = max(lims[3], arg.range[3]), min(lims[4], arg.range[4])
	row_dim, col_dim = (arg.layout == "" || arg.layout[2] == 'C') ? (1,2) : (2,1)	# If RowMajor the array is transposed
	slope = (size(arg, col_dim) - 1) / (arg.x[end] - arg.x[1])
	pix_x = round.(Int, slope .* (lims[1:2] .- arg.x[1]) .+ 1)
	slope = (size(arg, row_dim) - 1) / (arg.y[end] - arg.y[1])
	pix_y = round.(Int, slope .* (lims[3:4] .- arg.y[1]) .+ 1)

	function rearrange_ranges(pix_x, pix_y)
		# Rearrange the cropping limits if the layout is Rowmajor and/or Topdown
		if (arg.layout[1] == 'T')  pix_y = [size(arg, row_dim)-pix_y[2]+1, size(arg, row_dim)-pix_y[1]+1]	end
		if (arg.layout[2] == 'R')  pix_x, pix_y = pix_y, pix_x  end
		pix_x, pix_y
	end

	x, y = arg.x[pix_x[1]:pix_x[2]+1], arg.y[pix_y[1]:pix_y[2]+1]
	if (arg.layout != "")  pix_x, pix_y = rearrange_ranges(pix_x, pix_y)  end
	cropped = (ndims(arg) == 2) ? arg[pix_y[1]:pix_y[2], pix_x[1]:pix_x[2]] : arg[pix_y[1]:pix_y[2], pix_x[1]:pix_x[2], :]
	range = copy(arg.range)
	range[1:4] = [x[1], x[end], y[1], y[end]]
	if (eltype(arg) <: AbstractFloat)
		zmin, zmax = extrema_nan(cropped)
		range[5:6] = [zmin, zmax]
	end
	out = isa(arg, GMTgrid) ?  mat2grid(cropped, arg) : mat2img(cropped, arg)
	out.x, out.y, out.range = x, y, range
	out, pix_x, pix_y
end


text:
# This file is a part of LegendGeSim.jl, licensed under the MIT License (MIT).

using LegendGeSim
using Test

using LinearAlgebra, Random
# using Some, Other, Packages

@testset "utils" begin
    @testset "dummy_function" begin
        @test LegendGeSim.dummy_function(-4) == 16
        @test LegendGeSim.dummy_function(4) == 16
    end
end


text:
# SPDX-License-Identifier: MIT

const x = "//\""

### line1

const y = "/**\""
const c = 'c'

const v = v"v1.1.1"

const b = b"010001"

const raw = raw"xxx"

#=
 # line1
 # line2
 # line3
 =#

max_stars_repo_name: lxmwq/apidoc


text:
using PolygonMinkowskiSum
using Test

using GeometryTypes: Point
using PolygonMinkowskiSum: SimplePolygon, cross, find_ccw_violations

sweepme = SimplePolygon{Array{Point{2,Float64},1}}(Point{2,Float64}[[-0.35, -0.35], [0.35, -0.35], [0.35, 0.35], [-0.35, 0.35]])
alongme = SimplePolygon{Array{Point{2,Float64},1}}(Point{2,Float64}[[0.19857737682165164, 0.0], [0.19857737682165164, 0.14427490936067594], [0.06452170095561843, 0.19857737682165164], [-0.06452170095561849, 0.19857737682165164], [-0.19857737682165166, 0.1442749093606759], [-0.19857737682165166, 0.0], [-0.19857737682165166, -0.14427490936067594], [-0.06452170095561849, -0.19857737682165166], [0.06452170095561843, -0.19857737682165166], [0.19857737682165164, -0.144274909360676]])

@test find_ccw_violations(sweepme) == []
@test find_ccw_violations(alongme) == []

result = minkowski_sum(sweepme, alongme)

@test find_ccw_violations(result) == []


text:
@testset "vernal point" begin
    hours = 0:24
    lat,long = 0.0,0.0
    deg2second = 24*3600/360
    datetimes = DateTime(2021,3,20) .+ Hour.(hours) .- Second(round(long*deg2second))
    res2 = @pipe calc_sun_position_MOD.(datetime2julian.(datetimes)) |> DataFrame(_)
    # latitude over ecliptic of sun is zero
    @test all(res2.β .≈ 0.0)
    au2m = 149597870700.0
    # distance to sun about one astromical unit (150_000 km)
    @test all(.≈(res2.r,  au2m, rtol = 0.02))
    # declination crossing from negative to positive
    @test res2.δ[1] < 0
    @test res2.δ[end] > 0
    # azimuth crossing zero longitude
    @test res2.α[1] < 0
    @test res2.α[end] > 0
    # function tmpf()
    #     using Plots, StatsPlots
    #     @df res2 scatter(hours, cols([:α, :δ]), legend = :topleft, xlab="hours")
    #     @df res3 scatter(hours, cols(1:2), legend = :topleft, xlab="hours")
    # end
    #
    res3 = @pipe calc_sun_position_hor.(datetimes, lat, long) |> DataFrame(_)
    # altitude (-pi, +pi)
    @test minimum(res3.altitude) ≈ -π/2 atol = 0.1
    @test maximum(res3.altitude) ≈ +π/2 atol = 0.1
    #
    lat,long = 45.0,0.0
    datetimes = DateTime(2021,3,20) .+ Hour.(hours) .- Second(round(long*deg2second))
    res3 = @pipe calc_sun_position_hor.(datetimes, lat, long) |> DataFrame(_)
    # altitude (-pi, +pi)
    @test minimum(res3.altitude) ≈ -π/4 atol = 0.1
    @test maximum(res3.altitude) ≈ +π/4 atol = 0.1
    # azimuth increasing (aside edge cases)
    @test all(diff(res3.azimuth[2:(end-1)]) .> 0)
end

@testset "Dresden summer" begin
    hours = 0:24
    lat,long = 51.0, 13.6
    deg2second = 24*3600/360
    doy = 160
    datetimes = DateTime(2021) .+Day(doy-1) .+ Hour.(hours) .- Second(round(long*deg2second))
    # function tmpf()
    #     using Plots, StatsPlots
    #     @df res2 scatter(hours, cols([:α, :δ]), legend = :topleft, xlab="hours")
    #     @df res3 scatter(hours, cols(1:2), legend = :topleft, xlab="hours")
    # end
    #
    res3 = @pipe calc_sun_position_hor.(datetimes, lat, long) |> DataFrame(_)
    # maximum altitude at noon
    @test argmax(res3.altitude) == 1 + (length(hours)-1) ÷ 2
    # azimuth increasing (aside edge cases)
    @test all(diff(res3.azimuth[2:(end-1)]) .> 0)
end


text:
#module map_covid19
using Dates, DataFrames, ColorSchemes, PlotlyJS

#include("COVID19.jl")
#include("load_covid19_data.jl")

reasonable_resolution() = (1000, 800)

# time index
#for ti = 8:length(covid19us[1].confirmed)
#ti = 200
ti = length(covid19us[1].confirmed)

cconfirmed = [covid19us[i].confirmed[ti] for i in 1:length(covid19us)];
dcconfirmed = [mean(covid19us[i].confirmed[ti-7+1:ti] - covid19us[i].confirmed[ti-8+1:ti-1]) for i in 1:length(covid19us)];
dcconfirmedinst = [mean(covid19us[i].confirmed[ti] - covid19us[i].confirmed[ti-1]) for i in 1:length(covid19us)];
cdeath = [covid19us[i].death[ti] for i in 1:length(covid19us)];
dcdeath = [mean(covid19us[i].death[ti-7+1:ti] - covid19us[i].death[ti-8+1:ti-1]) for i in 1:length(covid19us)];
dcdeathinst = [mean(covid19us[i].death[ti] - covid19us[i].death[ti-1]) for i in 1:length(covid19us)];

# extracting county indices with land area and population
cind = findall((countyarea .> 0) .& (cpop .> 0));

tind = 1:length(covid19us[1].confirmed);
t = covid19us[1].time[tind];

# plot a map of confirmed cases
#cval = cconfirmed[cind] ./ (countyarea[cind] .* cpop[cind]);
cval = cconfirmed[cind] ./ cpop[cind];
log10cval = log10.(cval);
#sval = countyarea[cind] .* cpop[cind];
sval = cpop[cind];
log10sval = log10.(sval);

# daily confirmed cases per 100k averaged over 7 days
dCONFIRMEDpc = dcconfirmed[cind] ./ cpop[cind] .* 1e5;
bind = findall(dCONFIRMEDpc .<= 0);
dCONFIRMEDpc[bind] .= NaN;

# daily instant confirmed cases per 100k
dCONFIRMEDpcinst = dcconfirmedinst[cind] ./ cpop[cind] .* 1e5;
bind = findall(dCONFIRMEDpcinst .<= 0);
dCONFIRMEDpcinst[bind] .= NaN;

# daily death per 100k averaged over 7 days
dDEATHpc = dcdeath[cind] ./ cpop[cind] .* 1e5;
bind = findall(dDEATHpc .<= 0);
dDEATHpc[bind] .= NaN;

# daily instant death per 100k
dDEATHpcinst = dcdeathinst[cind] ./ cpop[cind] .* 1e5;
bind = findall(dDEATHpcinst .<= 0);
dDEATHpcinst[bind] .= NaN;

df = DataFrame(LON = clon[cind],
    LAT = clat[cind],
    POPULATION = cpop[cind],
    log10POPULATION = log10.(cpop[cind]),
    logPOPULATION = log.(cpop[cind]),
    nrootPOPULATION = (cpop[cind]) .^ (1.0/2.85)/5,

    CONFIRMEDpc = cconfirmed[cind] ./ cpop[cind] .* 1e5,
    CONFIRMED = cconfirmed[cind],
    log10CONFIRMEDpc = log10.(cconfirmed[cind] ./ cpop[cind] .* 1e5),

    dCONFIRMEDpc = dCONFIRMEDpc,
    dCONFIRMED = dcconfirmed[cind],
    log10dCONFIRMEDpc = log10.(dCONFIRMEDpc),

    dCONFIRMEDpcinst = dCONFIRMEDpcinst,
    dCONFIRMEDinst = dcconfirmedinst[cind],
    log10dCONFIRMEDpcinst = log10.(dCONFIRMEDpcinst),

    DEATHpc = cdeath[cind] ./ cpop[cind] .* 1e5,
    DEATH = cdeath[cind],
    log10DEATHpc = log10.(cdeath[cind] ./ cpop[cind] .* 1e5),

    dDEATHpc = dDEATHpc,
    dDEATH = dcdeath[cind],
    log10dDEATHpc = log10.(dDEATHpc),

    dDEATHpcinst = dDEATHpcinst,
    dDEATHinst = dcdeathinst[cind],
    log10dDEATHpcinst = log10.(dDEATHpcinst),
    );

function pmapconfirmed()
    pmapsg = PlotlyJS.scattergeo(; locationmode="USA-states",
        lat = df[!, :LAT],
        lon = df[!, :LON],
        marker_size=df[!, :nrootPOPULATION],
        marker_color=df[!, :log10CONFIRMEDpc],
        marker_colorscale="Jet",
        marker_showscale = true,
        marker_cmin = 3,
        marker_cmax = 4.3,

        hoverinfo="text",
        hovertext=
            "County: " .* county[cind] .* ", " .* state[cind] .* "
" .*
            "Population: " .* string.(df[!, :POPULATION]) .* "
" .*
            "Confirmed Cases (actual): " .* string.(df[!, :CONFIRMED]) .* "
" .*
            "Confirmed Cases per 100k: " .* string.(round.(df[!, :CONFIRMEDpc], digits=2)) .* "
" .*
            "Cumulative Deaths (actual): " .* string.(df[!, :DEATH]) .* "
" .*
            "Cumulative Deaths per 100k: " .* string.(round.(df[!, :DEATHpc], digits=2)) .* "
" .*
            "% infected with COVID19: " .* string.(round.(df[!, :CONFIRMED] ./ df[!, :POPULATION] .* 100, digits=2)) .* "
" .*
            "% death rate: " .* string.(round.(df[!, :DEATH] ./ df[!, :CONFIRMED] .*100, digits=2))
        );

    geo = attr(scope="usa",
        projection_type="ablbers usa",
        showland=true,
        landcolor="rgb(220, 220, 220)",
        subunitwidth=1,
        countrywidth=1,
        subunitcolor="rgb(255,255,255)",
        countrycolor="rgb(255,255,255)",
        );

    title = attr(text = "Cumulative Confirmed Cases per 100k " * string(t[ti])[1:10] * " (log10)
https://truedichotomy.github.io/covid19_public/
", yref = "paper", y=1.0)

    layout = Layout(; title=title, showlegend=false, geo=geo)
    PlotlyJS.plot(pmapsg, layout)
end

function pmapdconfirmed()
    pmapsg = PlotlyJS.scattergeo(; locationmode="USA-states",
        lat = df[!, :LAT],
        lon = df[!, :LON],
        #marker_size=df[!, :logPOPULATION]/1.3,
        marker_size=df[!, :nrootPOPULATION],
        marker_color=df[!, :log10dCONFIRMEDpc],
        marker_colorscale="Jet",
        marker_showscale = true,
        marker_cmin = -0.5,
        marker_cmax = 2.5,
        hoverinfo="text",
        #hovertext=string.(df[!, :dCONFIRMED]) .* " " .* county[cind],
        hovertext=
            "County: " .* county[cind] .* ", " .* state[cind] .* "
" .*
            "Population: " .* string.(df[!, :POPULATION]) .* "
" .*
            "Confirmed Cases (actual): " .* string.(df[!, :CONFIRMED]) .* "
" .*
            "Confirmed Cases per 100k: " .* string.(round.(df[!, :CONFIRMEDpc], digits=2)) .* "
" .*
            "Daily New Cases (actual): " .* string.(round.(df[!, :dCONFIRMED], digits=2)) .* "
" .*
            "Daily New Cases per 100k: " .* string.(round.(df[!, :dCONFIRMEDpc], digits=2)) .* "
" .*
            "% infected with COVID19: " .* string.(round.(df[!, :CONFIRMED] ./ df[!, :POPULATION] .* 100, digits=2)) .* "
" .*
            "% death rate: " .* string.(round.(df[!, :DEATH] ./ df[!, :CONFIRMED] .*100, digits=2))
        );

    geo = attr(scope="usa",
        projection_type="ablbers usa",
        showland=true,
        landcolor="rgb(220, 220, 220)",
        subunitwidth=1,
        countrywidth=1,
        subunitcolor="rgb(255,255,255)",
        countrycolor="rgb(255,255,255)",
        );

    title = attr(text = "Daily New Cases per 100k " * string(t[ti])[1:10] * " (log10, 7 day avg.)
https://truedichotomy.github.io/covid19_public/
", yref = "paper", y=1.0)

    layout = Layout(; title=title, showlegend=false, geo=geo)
    PlotlyJS.plot(pmapsg, layout)
end

function pmapdconfirmedinst()
    pmapsg = PlotlyJS.scattergeo(; locationmode="USA-states",
        lat = df[!, :LAT],
        lon = df[!, :LON],
        marker_size=df[!, :nrootPOPULATION],
        #marker_size=df[!, :logPOPULATION]/1.3,
        marker_color=df[!, :log10dCONFIRMEDpcinst],
        marker_colorscale="Jet",
        marker_showscale = true,
        marker_cmin = -0.5,
        marker_cmax = 2.5,
        hoverinfo="text",
        #hovertext=string.(df[!, :dCONFIRMED]) .* " " .* county[cind],
        hovertext=
            "County: " .* county[cind] .* ", " .* state[cind] .* "
" .*
            "Population: " .* string.(df[!, :POPULATION]) .* "
" .*
            "Confirmed Cases (actual): " .* string.(df[!, :CONFIRMED]) .* "
" .*
            "Confirmed Cases per 100k: " .* string.(round.(df[!, :CONFIRMEDpc], digits=2)) .* "
" .*
            "Daily New Cases (actual): " .* string.(df[!, :dCONFIRMEDinst]) .* "
" .*
            "Daily New Cases per 100k: " .* string.(round.(df[!, :dCONFIRMEDpcinst], digits=2)) .* "
" .*
            "% infected with COVID19: " .* string.(round.(df[!, :CONFIRMED] ./ df[!, :POPULATION] .* 100, digits=2)) .* "
" .*
            "% death rate: " .* string.(round.(df[!, :DEATH] ./ df[!, :CONFIRMED] .*100, digits=2))
        );

    geo = attr(scope="usa",
        projection_type="ablbers usa",
        showland=true,
        landcolor="rgb(220, 220, 220)",
        subunitwidth=1,
        countrywidth=1,
        subunitcolor="rgb(255,255,255)",
        countrycolor="rgb(255,255,255)",
        );

    title = attr(text = "Daily New Cases per 100k " * string(t[ti])[1:10] * " (log10)
https://truedichotomy.github.io/covid19_public/
", yref = "paper", y=1.0)

    layout = Layout(; title=title, showlegend=false, geo=geo)
    PlotlyJS.plot(pmapsg, layout)
end

function pmapdeath()
    pmapsg = PlotlyJS.scattergeo(; locationmode="USA-states",
        lat = df[!, :LAT],
        lon = df[!, :LON],
        marker_size=df[!, :nrootPOPULATION],
        #marker_size=df[!, :logPOPULATION]/1.3,
        marker_color=df[!, :log10DEATHpc],
        marker_colorscale="Jet",
        marker_showscale = true,
        marker_cmin = 0,
        marker_cmax = 2.5,
        hoverinfo="text",
        #hovertext=string.(df[!, :DEATH]) .* " " .* county[cind],
        hovertext=
            "County: " .* county[cind] .* ", " .* state[cind] .* "
" .*
            "Population: " .* string.(df[!, :POPULATION]) .* "
" .*
            "Confirmed Cases (actual): " .* string.(df[!, :CONFIRMED]) .* "
" .*
            "Confirmed Cases per 100k: " .* string.(round.(df[!, :CONFIRMEDpc], digits=2)) .* "
" .*
            "Cumulative Deaths (actual): " .* string.(df[!, :DEATH]) .* "
" .*
            "Cumulative Deaths per 100k: " .* string.(round.(df[!, :DEATHpc], digits=2)) .* "
" .*
            "% infected with COVID19: " .* string.(round.(df[!, :CONFIRMED] ./ df[!, :POPULATION] .* 100, digits=2)) .* "
" .*
            "% death rate: " .* string.(round.(df[!, :DEATH] ./ df[!, :CONFIRMED] .*100, digits=2))
        );

    geo = attr(scope="usa",
        projection_type="ablbers usa",
        showland=true,
        landcolor="rgb(220, 220, 220)",
        subunitwidth=1,
        countrywidth=1,
        subunitcolor="rgb(255,255,255)",
        countrycolor="rgb(255,255,255)",
        );

    title = attr(text = "Cumulative Death per 100k " * string(t[ti])[1:10] * " (log10)
https://truedichotomy.github.io/covid19_public/
", yref = "paper", y=1.0)

    layout = Layout(; title=title, showlegend=false, geo=geo)
    PlotlyJS.plot(pmapsg, layout)
end

function pmapddeath()
    pmapsg = PlotlyJS.scattergeo(; locationmode="USA-states",
        lat = df[!, :LAT],
        lon = df[!, :LON],
        #marker_size=df[!, :logPOPULATION]/1.3,
        marker_size=df[!, :nrootPOPULATION],
        marker_color=df[!, :log10dDEATHpc],
        marker_colorscale="Jet",
        marker_showscale = true,
        marker_cmin = -1.5,
        marker_cmax = 1,
        hoverinfo="text",
        #hovertext=string.(df[!, :dDEATH]) .* " " .* county[cind],
        hovertext=
            "County: " .* county[cind] .* ", " .* state[cind] .* "
" .*
            "Population: " .* string.(df[!, :POPULATION]) .* "
" .*
            "Total Deaths (actual): " .* string.(df[!, :DEATH]) .* "
" .*
            "Total Deaths per 100k: " .* string.(round.(df[!, :DEATHpc], digits=2)) .* "
" .*
            "Daily New Deaths (actual): " .* string.(round.(df[!, :dDEATH], digits=2)) .* "
" .*
            "Daily New Deaths per 100k: " .* string.(round.(df[!, :dDEATHpc], digits=2)) .* "
" .*
            "% infected with COVID19: " .* string.(round.(df[!, :CONFIRMED] ./ df[!, :POPULATION] .* 100, digits=2)) .* "
" .*
            "% death rate: " .* string.(round.(df[!, :DEATH] ./ df[!, :CONFIRMED] .*100, digits=2))
        );

    geo = attr(scope="usa",
        projection_type="ablbers usa",
        showland=true,
        landcolor="rgb(220, 220, 220)",
        subunitwidth=1,
        countrywidth=1,
        subunitcolor="rgb(255,255,255)",
        countrycolor="rgb(255,255,255)",
        );

    title = attr(text = "Daily New Deaths per 100k " * string(t[ti])[1:10] * " (log10, 7 day avg.)
https://truedichotomy.github.io/covid19_public/
", yref = "paper", y=1.0)

    layout = Layout(; title=title, showlegend=false, geo=geo)
    PlotlyJS.plot(pmapsg, layout)
end

function pmapddeathinst()
    pmapsg = PlotlyJS.scattergeo(; locationmode="USA-states",
        lat = df[!, :LAT],
        lon = df[!, :LON],
        #marker_size=df[!, :logPOPULATION]/1.3,
        marker_size=df[!, :nrootPOPULATION],

        marker_color=df[!, :log10dDEATHpcinst],
        marker_colorscale="Jet",
        marker_showscale = true,
        marker_cmin = -1.5,
        marker_cmax = 1,
        hoverinfo="text",
        #hovertext=string.(df[!, :dDEATH]) .* " " .* county[cind],
        hovertext=
            "County: " .* county[cind] .* ", " .* state[cind] .* "
" .*
            "Population: " .* string.(df[!, :POPULATION]) .* "
" .*
            "Total Deaths (actual): " .* string.(df[!, :DEATH]) .* "
" .*
            "Total Deaths per 100k: " .* string.(round.(df[!, :DEATHpc], digits=2)) .* "
" .*
            "Daily New Deaths (actual): " .* string.(df[!, :dDEATHinst]) .* "
" .*
            "Daily New Deaths per 100k: " .* string.(round.(df[!, :dDEATHpcinst], digits=2)) .* "
" .*
            "% infected with COVID19: " .* string.(round.(df[!, :CONFIRMED] ./ df[!, :POPULATION] .* 100, digits=2)) .* "
" .*
            "% death rate: " .* string.(round.(df[!, :DEATH] ./ df[!, :CONFIRMED] .*100, digits=2))
        );

    geo = attr(scope="usa",
        projection_type="ablbers usa",
        showland=true,
        landcolor="rgb(220, 220, 220)",
        subunitwidth=1,
        countrywidth=1,
        subunitcolor="rgb(255,255,255)",
        countrycolor="rgb(255,255,255)",
        );

    title = attr(text = "Daily New Deaths per 100k " * string(t[ti])[1:10] * " (log10)
https://truedichotomy.github.io/covid19_public/
", yref = "paper", y=1.0)

    layout = Layout(; title=title, showlegend=false, geo=geo)
    PlotlyJS.plot(pmapsg, layout)
end

#figoutdir = "/Volumes/GoogleDrive/My Drive/COVID19/";
figoutdir = "../covid19_public/maps/";

#PlotlyJS._js_path = "https://cdn.plot.ly/plotly-latest.min.js";

PlotlyJS.savefig(pmapconfirmed(), figoutdir * "covid19map_confirmed.html", format="html");
PlotlyJS.savefig(pmapdconfirmed(), figoutdir * "covid19map_delta_confirmed_7days.html", format="html");
PlotlyJS.savefig(pmapdconfirmedinst(), figoutdir * "covid19map_delta_confirmed_latest.html", format="html");
PlotlyJS.savefig(pmapdeath(), figoutdir * "covid19map_death.html", format="html");
PlotlyJS.savefig(pmapddeath(), figoutdir * "covid19map_delta_death_7days.html", format="html");
PlotlyJS.savefig(pmapddeathinst(), figoutdir * "covid19map_delta_death_latest.html", format="html");

#PlotlyJS.savehtml(pmapconfirmed(), figoutdir * "covid19map_confirmed.html", :remote);
#PlotlyJS.savehtml(pmapdconfirmed(), figoutdir * "covid19map_delta_confirmed_7days.html", :remote);
#PlotlyJS.savehtml(pmapdconfirmedinst(), figoutdir * "covid19map_delta_confirmed_latest.html", :remote);
#PlotlyJS.savehtml(pmapdeath(), figoutdir * "covid19map_death.html", :remote);
#PlotlyJS.savehtml(pmapddeath(), figoutdir * "covid19map_delta_death_7days.html", :remote);
#PlotlyJS.savehtml(pmapddeathinst(), figoutdir * "covid19map_delta_death_latest.html", :remote);
#end


text:
using Test
using ConvexFit

using ConvexFit: SolverState, SolverTrace, update!, Cache, CVCache, LOOCVCache
using Optim: Brent, optimize

# Threads.@threads may cause the CI to freeze on Windows in older versions
const testthreads = VERSION >= v"1.2" && !Sys.iswindows() || VERSION >= v"1.7"

function fmin(f::Function)
    r = optimize(f, 0.0, 100.0, Brent(), abs_tol=1e-4, store_trace=true)
    return r, r.minimizer
end

const tests = [
    "fit",
    "select"
]

printstyled("Running tests:\n", color=:blue, bold=true)

@time for test in tests
    include("$test.jl")
    println("\033[1m\033[32mPASSED\033[0m: $(test)")
end


text:
import Base.@deprecate, Base.@deprecate_binding
import Base.depwarn

# DEPRECATED: defined for compatibility
type NucleotideSequence{T<:Nucleotide} end

@deprecate NucleotideSequence BioSequence
@deprecate NucleotideSequence(::Type{DNANucleotide}) DNASequence
@deprecate NucleotideSequence(::Type{RNANucleotide}) RNASequence
@compat begin
    (::Type{NucleotideSequence{DNANucleotide}})() = DNASequence()
    (::Type{NucleotideSequence{RNANucleotide}})() = RNASequence()
    (::Type{NucleotideSequence{DNANucleotide}})(
        seq::Union{AbstractVector{DNANucleotide},AbstractString}) = DNASequence(seq)
    (::Type{NucleotideSequence{RNANucleotide}})(
        seq::Union{AbstractVector{RNANucleotide},AbstractString}) = RNASequence(seq)
end
Base.convert(::Type{NucleotideSequence}, seq::DNAKmer) = DNASequence(seq)
Base.convert(::Type{NucleotideSequence}, seq::RNAKmer) = RNASequence(seq)
Base.convert(::Type{NucleotideSequence{DNANucleotide}}, seq::Union{AbstractVector{DNANucleotide},AbstractString,DNAKmer}) = DNASequence(seq)
Base.convert(::Type{NucleotideSequence{RNANucleotide}}, seq::Union{AbstractVector{RNANucleotide},AbstractString,RNAKmer}) = RNASequence(seq)

@deprecate NucleotideCounts(seq) Composition(seq)
@deprecate npositions(seq::BioSequence) ambiguous_positions(seq)

@deprecate kmer Kmer
@deprecate dnakmer DNAKmer
@deprecate rnakmer RNAKmer


text:
using DataFrames
using Feather
using Distributions

# params
LAMBDA = 0.6882627
ALPHA = 1.3441739239406807
SIGMA = 0.4203008

function bias_lambda(slice, alpha=ALPHA, sd=SIGMA)
    sd * (- alpha * sin(2 * pi  * slice))
end

function draw_arrival_rate(slice, n, arrival_rate_mean=LAMBDA; alpha=ALPHA, arrival_rate_sd=SIGMA)
    rel_slice = slice / n
    lambda = arrival_rate_mean + bias_lambda(rel_slice, alpha, arrival_rate_sd)
    dist = Poisson(lambda)
    return rand(dist)
end

# one repetition simulation
n = 1440
average_day = [(i, draw_arrival_rate(i, n)) for i in 0:n-1]
average_day_df = DataFrame(average_day)
rename!(average_day_df, [1 => :time_index, 2 => :arr_count])
average_day_df[:, :hour_of_day] = [i for i in 0:23 for j in 1:60]
Feather.write(joinpath("data", "simulated_arrivals.feather"), average_day_df)

# bootstrap simulation (100 repititions)
n = 1440
average_day = [(i, mean([draw_arrival_rate(i, n) for j in 1:100])) for i in 0:n-1 ]
average_day_df = DataFrame(average_day)
rename!(average_day_df, [1 => :time_index, 2 => :arr_count])
average_day_df[:, :hour_of_day] = [i for i in 0:23 for j in 1:60]
Feather.write(joinpath("data", "simulated_arrivals.feather"), average_day_df)

max_stars_repo_name: JohannesNakayama/hacker-news-analysis

max_issues_repo_name: JohannesNakayama/hacker-news-analysis

max_forks_repo_name: JohannesNakayama/hacker-news-analysis


text:
macro stack(expr...)
	esc(genstack(collect(expr)))
end

function genstack(e::Vector)
	q  = Expr(:new, :List , nothing  )
	for i in e
		q = Expr(:new, :List , i , q  )
	end
	return q
end

max_stars_repo_name: frankill/ftest

max_issues_repo_name: frankill/ftest

max_forks_repo_name: frankill/ftest


text:
using FractionalCalculus
using Plots, LaTeXStrings

target = gamma(6)/gamma(2.4)*tspan.^1.4

s="\$D^{3.6}x^5\$"

tspan=collect(0:0.01:6)

result=fracdiff(x->x^5, 3.6, 6, 0.01, RLDiffMatrix())

plot(tspan, result, title=s, legend=:bottomright, label="Numerical")
plot!(tspan, target, lw=3, ls=:dash, label="Analytical")
savefig("./arbitrary_order_derivative.png")


text:
using CUDex

inputsize = 16
seqlength = 10
batchsize = 4
hiddensize = 32
numlayers = 1
dropout = 0.5

x = DexArray(rand(seqlength,batchsize,inputsize)) #

y,hy,cy,w,rd,rs = rnn_ft(x,hiddensize=16)
display(y)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

@enum CONTEXT_TYPE CPU=1 GPU=2 CPU_PINNED=3

Base.convert(::Type{CONTEXT_TYPE}, x::Integer) = CONTEXT_TYPE(x)

"""
    Context(dev_type, dev_id)

A context describes the device type and id on which computation should be carried on.
"""
struct Context
  device_type::CONTEXT_TYPE
  device_id::Int

  Context(dev_type::CONTEXT_TYPE, dev_id::Integer = 0) = new(dev_type, dev_id)
end

Context(dev_type::Integer, dev_id::Integer = 0) =
  Context(convert(CONTEXT_TYPE, dev_type), dev_id)

Base.show(io::IO, ctx::Context) =
  print(io, "$(ctx.device_type)$(ctx.device_id)")

"""
    cpu(dev_id)

Get a CPU context with a specific id. `cpu()` is usually the default context for many
operations when no context is specified.

# Arguments
* `dev_id::Integer = 0`: the CPU id.
"""
cpu(dev_id::Integer = 0) = Context(CPU, dev_id)

"""
    gpu(dev_id)

Get a GPU context with a specific id. The K GPUs on a node is typically numbered as 0,...,K-1.

# Arguments
* `dev_id::Integer = 0` the GPU device id.
"""
gpu(dev_id::Integer = 0) = Context(GPU, dev_id)

"""
    num_gpus()

Query CUDA for the number of GPUs present.
"""
function num_gpus()
  n = Ref{Cint}()
  @mxcall :MXGetGPUCount (Ref{Cint},) n
  n[]
end

"""
    gpu_memory_info(dev_id = 0)::Tuple{UInt64,UInt64}

Query CUDA for the free and total bytes of GPU global memory.
It returns a tuple of `(free memory, total memory)`.

```julia-repl
julia> mx.gpu_memory_info()
(0x00000003af240000, 0x00000003f9440000)
```
"""
function gpu_memory_info(dev_id = 0)
  free = Ref{UInt64}()
  n = Ref{UInt64}()
  @mxcall :MXGetGPUMemoryInformation64 (Cint, Ref{UInt64}, Ref{UInt64}) dev_id free n
  free[], n[]
end


text:
using SafeTestsets

module GuidingCenter4dTests

    using Test
    using GeometricIntegrators

    const Δt = 1.
    const nt = 1

    const nl = 100
    const nx = 10
    const ny = 10

    export test_guiding_center_4d_glrk, test_guiding_center_4d_vpglrk
    export nl, nx, ny

    function test_guiding_center_4d_glrk(ode; Δt = Δt)
        int = Integrator(ode, getTableauGLRK(1), Δt)
        sol = integrate(ode, int, nt)
        @test true
    end

    function test_guiding_center_4d_vpglrk(ode; Δt = Δt)
        int = Integrator(ode, getTableauVPGLRK(1), Δt)
        sol = integrate(ode, int, nt)
        @test true
    end

end

@safetestset "Guiding Centre Dynamics in 4D with ITER-like Solov'ev Equilibrium with X-Point                      " begin

    using ChargedParticleDynamics.GuidingCenter4d.SolovevIterXpoint
    using ..GuidingCenter4dTests

    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_trapped()...), Δt=10000.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_barely_passing()...), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_barely_trapped()...), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_deeply_passing()...), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_deeply_trapped()...), Δt=1.)
    # test_guiding_center_4d_glrk(guiding_center_4d_loop_ode(nl), Δt=1.)
    # test_guiding_center_4d_glrk(guiding_center_4d_surface_ode(nx, ny), Δt=1.)

    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_trapped()...), Δt=10000.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_barely_passing()...), Δt=1.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_barely_trapped()...), Δt=1.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_deeply_passing()...), Δt=1.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_deeply_trapped()...), Δt=1.)
    # test_guiding_center_4d_vpglrk(guiding_center_4d_loop_iode(nl), Δt=1.)
    # test_guiding_center_4d_vpglrk(guiding_center_4d_surface_iode(nx, ny), Δt=1.)

end

@safetestset "Guiding Centre Dynamics in 4D with medium-size Tokamak Equilibrium in Cartesian Coordinates         " begin

    using ChargedParticleDynamics.GuidingCenter4d.TokamakMediumCartesian
    using ..GuidingCenter4dTests

    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_barely_passing()...), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_barely_trapped()...), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_deeply_passing()...), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_deeply_trapped()...), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_loop_ode(nl), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_surface_ode(nx, ny), Δt=1.)

    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_barely_passing()...), Δt=1.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_barely_trapped()...), Δt=1.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_deeply_passing()...), Δt=1.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_deeply_trapped()...), Δt=1.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_loop_iode(nl), Δt=.1)
    test_guiding_center_4d_vpglrk(guiding_center_4d_surface_iode(nx, ny), Δt=.1)

end

@safetestset "Guiding Centre Dynamics in 4D with medium-size Tokamak Equilibrium in Cylindrical Coordinates       " begin

    using ChargedParticleDynamics.GuidingCenter4d.TokamakMediumCylindrical
    using ..GuidingCenter4dTests

    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_barely_passing()...), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_barely_trapped()...), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_deeply_passing()...), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_deeply_trapped()...), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_loop_ode(nl), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_surface_ode(nx, ny), Δt=1.)

    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_barely_passing()...), Δt=1.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_barely_trapped()...), Δt=1.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_deeply_passing()...), Δt=1.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_deeply_trapped()...), Δt=1.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_loop_iode(nl), Δt=.1)
    test_guiding_center_4d_vpglrk(guiding_center_4d_surface_iode(nx, ny), Δt=.1)

end

@safetestset "Guiding Centre Dynamics in 4D with small-size Tokamak Equilibrium in Cartesian Coordinates          " begin

    using ChargedParticleDynamics.GuidingCenter4d.TokamakSmallCartesian
    using ..GuidingCenter4dTests

    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_barely_passing()...), Δt=400.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_barely_trapped()...), Δt=400.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_deeply_passing()...), Δt=400.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_deeply_trapped()...), Δt=400.)
    test_guiding_center_4d_glrk(guiding_center_4d_loop_ode(nl), Δt=400.)
    test_guiding_center_4d_glrk(guiding_center_4d_surface_ode(nx, ny), Δt=400.)

    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_barely_passing()...), Δt=10.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_barely_trapped()...), Δt=10.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_deeply_passing()...), Δt=10.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_deeply_trapped()...), Δt=10.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_loop_iode(nl), Δt=10.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_surface_iode(nx, ny), Δt=10.)

end

@safetestset "Guiding Centre Dynamics in 4D with small-size Tokamak Equilibrium in Cylindrical Coordinates        " begin

    using ChargedParticleDynamics.GuidingCenter4d.TokamakSmallCylindrical
    using ..GuidingCenter4dTests

    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_barely_passing()...), Δt=400.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_barely_trapped()...), Δt=400.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_deeply_passing()...), Δt=400.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_deeply_trapped()...), Δt=400.)
    test_guiding_center_4d_glrk(guiding_center_4d_loop_ode(nl), Δt=400.)
    test_guiding_center_4d_glrk(guiding_center_4d_surface_ode(nx, ny), Δt=400.)

    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_barely_passing()...), Δt=10.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_barely_trapped()...), Δt=10.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_deeply_passing()...), Δt=10.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_deeply_trapped()...), Δt=10.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_loop_iode(nl), Δt=10.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_surface_iode(nx, ny), Δt=10.)

end

@safetestset "Guiding Centre Dynamics in 4D with small-size Tokamak Equilibrium in Toroidal Coordinates           " begin

    using ChargedParticleDynamics.GuidingCenter4d.TokamakSmallToroidal
    using ..GuidingCenter4dTests

    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_barely_passing()...), Δt=400.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_barely_trapped()...), Δt=400.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_deeply_passing()...), Δt=400.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_deeply_trapped()...), Δt=400.)
    test_guiding_center_4d_glrk(guiding_center_4d_loop_ode(nl), Δt=400.)
    test_guiding_center_4d_glrk(guiding_center_4d_surface_ode(nx, ny), Δt=400.)

    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_barely_passing()...), Δt=400.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_barely_trapped()...), Δt=400.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_deeply_passing()...), Δt=400.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_deeply_trapped()...), Δt=400.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_loop_iode(nl), Δt=400.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_surface_iode(nx, ny), Δt=400.)

end

@safetestset "Guiding Centre Dynamics in 4D with symmetric Solov'ev Equilibrium                                   " begin

    using ChargedParticleDynamics.GuidingCenter4d.SolovevSymmetricField
    using ..GuidingCenter4dTests

    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_barely_passing()...), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_barely_trapped()...), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_deeply_passing()...), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_ode(initial_conditions_deeply_trapped()...), Δt=1.)

    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_barely_passing()...), Δt=1.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_barely_trapped()...), Δt=1.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_deeply_passing()...), Δt=1.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_iode(initial_conditions_deeply_trapped()...), Δt=1.)

end

@safetestset "Guiding Centre Dynamics in 4D with symmetric Equlibrium                                             " begin

    using ChargedParticleDynamics.GuidingCenter4d.SymmetricField
    using ..GuidingCenter4dTests

    test_guiding_center_4d_glrk(guiding_center_4d_loop_ode(nl), Δt=1.)
    test_guiding_center_4d_glrk(guiding_center_4d_surface_ode(nx, ny), Δt=1.)

    test_guiding_center_4d_vpglrk(guiding_center_4d_loop_iode(nl), Δt=1.)
    test_guiding_center_4d_vpglrk(guiding_center_4d_surface_iode(nx, ny), Δt=1.)

end

@safetestset "Guiding Centre Dynamics in 4D with Theta Pinch Equilibrium                                          " begin

    using ChargedParticleDynamics.GuidingCenter4d.ThetaPinchField
    using ..GuidingCenter4dTests

    test_guiding_center_4d_glrk(guiding_center_4d_loop_ode(nl), Δt=1.)

    test_guiding_center_4d_vpglrk(guiding_center_4d_loop_iode(nl), Δt=1.)

end


text:
#1. Important Includes:
#-------------------------------------------------------------------------------
println(workers())

@everywhere using SegwayRobot3# Include Environment
@everywhere using Knet
@everywhere using AutoGrad
@everywhere using PyCall

#using ModelVisualisations #To have videos
using ProgressMeter # To estimate the simulation time
using PyPlot        # To visualize the improvement
using JLD2 #To save the results

include("data_types.jl") #Required Data types
include("trajectory_generating.jl") #Trajectory and Episodes
include("net_functionality.jl") #Neuronal Network
include("parallel_execution.jl") #To run things in parallel
include("utils.jl") #Algo-Specifics
include("TRPO.jl") #Including 2nd order derivatives, line-seach and trust region
include("GAE.jl") #Include General Advantage Estimation
include("PPO.jl") #Include the Surrogate Objective calculation
include("visualization.jl") #Visualization
norm_active = true #set TRUE if normalization should be activated
include("layer_norm.jl") #Normalization

include("NIPS2018.jl") #If trained on the nips environment

#2: Initialize
#-------------------------------------------------------------------------------

#2.1: Environment Parameters
#------------------------
params           = Params()
params.dt        = 1/100;   #Timestep
params.t_horizon = 5;       #Time until end
params.l_traj    = length(0:params.dt:params.t_horizon); #Total timesteps per episode
params.u_range   = u_range    #Range of motor power (Scale of NN output)
params.s_init    = s_init; #Init_position: Pendulum on the bottom!
params.l_s_t     = length(params.s_init) #Contains number of relevant input states
params.l_out     = length(params.u_range) #Contains output dimension

#2.2: Parameters of the Network and Training
#------------------------
params.structure_net_pl  = Int[params.l_s_t + 0; 2^6; 2^6; params.l_out] #Net structure policy function
params.structure_net_vf  = Int[params.l_s_t + 0; 2^6; 2^6; 1] #Net structure value function
params.structure_net_qf  = Int[params.l_s_t + params.l_out; 2^8; 2^8; 1] #Net structure Q function 1
params.structure_net_qf2 = Int[params.l_s_t + params.l_out; 2^9; 2^7; 1] #Net structure Q function 2
params.init_range        = 0.001; #Range in which the weights are initialized
params.batch_size        = 16; #Size of number of timesteps for the Minibatch update
params.lr_pl             = 3e-4; #Learning rate for policy
params.lr_vf             = 3e-4; #Learning rate for value function
params.lr_qf             = 3e-4; #Learning rate for q function
params.σ                 = 0.05 * u_range; #Standard distribution for stochastic policy
params.n_traj            = 400; #Number of trajectories per iteration
params.n_iter            = 1500; #Number of training updates
params.buffer_size       = Int(floor(1e6 / params.batch_size)); #Maximum size of number of Minibatches in buffer
params.γ_disc            = 0.99; #Discount factor
params.λ_critic          = 0.0; #TD(λ) factor -0 is TD(0)-full history -1 is TD(1) only one step
params.λ_actor           = 0.95; #GAE factor
params.eps               = 0.2; #Clipping factor for PPO
params.δ                 = 0.01; #Delta Factor for TRPO update
params.cg_damping        = 0.1; #cg-damping factor for Fisher-Vector-Product
params.τ                 = 0.002; #Updating for running updates (mainly DDPG and SAC)
params.σ_off_policy      = 0.1; #SARSA regularization for off-policy q-training
params.clip_σ            = 0.5; #clipping for sarsa training
params.d_update          = 5; #Relation between q-function to policy training
#-------------------------------------------------------------------------------

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
# This file is pat of RLinearAlgebra.jl
# 1. Specifies type
# 2. Implements rsubsolve function
# 3. Exports Type

#suing LinearAlgebra

"""
    LinSysVecColProjFO <: LinSysVecColProjection

A mutable structure that represents the standard column projection method
with full orthogonalization against all previous projections.

See Patel, Vivak, Mohammad Jahangoshahi, and Daniel Adrian Maldonado. "An Implicit
    Representation and Iterative Solution of Randomly Sketched Linear Systems."
    SIAM Journal on Matrix Analysis and Applications (2021) 42:2, 800-831.
    https://doi.org/10.1137/19M1259481

# Fields

- `S::Union{Matrix{Float64}, Nothing}` is a matrix used for orthogonalizing against
    all previous search directions.

Calling `LinSysVecColProjFO()` defaults to `LinSysVecColProjFO(nothing)`.
"""
mutable struct LinSysVecColProjFO <: LinSysVecColProjection
    S::Union{Matrix{Float64}, Nothing}
end
LinSysVecColProjFO() = LinSysVecColProjFO(nothing)

# Common rsubsolve interface for linear systems
function rsubsolve!(
    type::LinSysVecColProjFO,
    x::AbstractVector,
    samp::Tuple{U,V,W} where {U<:AbstractVector,V<:AbstractArray,W<:Real},
    iter::Int64,
)
    # samp[1] is the search direction
    # samp[2] is the coefficient matrix A
    # samp[3] is the residual of normal system in A * samp[1], samp[1]' A' * (A * x - b)

    if iter == 1
        d = length(x)
        type.S = diagm(ones(Float64, d))
    end

    # Compute orthogonal component of A * samp[1] using S (projection matrix)
    q = samp[2]' * (samp[2] * samp[1])
    u = type.S * q

    if dot(u,u) < eps()^1.5 * length(x) #1.5 instead of 2 since normal system
        return nothing
    end

    # Update Iterate
    γ = dot(u, q)
    x .= x - u * (samp[3] / γ)

    # Update Projection Matrix
    type.S .= (I - ( u / γ) * q') * type.S

    return nothing
end


text:
module AdjacentFloat

!isdefined(Float) && typealias Float AbstractFloat;

export prevFloat, nextFloat, nFloatsSeparate, nextNearerToZero, nextAwayFromZero

include("adjacentFloats.jl")

end # module


text:
# Mixed Integer Problems
# from Kusum Deep, Krishna Pratap Singh, M. L. Kansal, and C. Mohan,
# "A real coded  genetic algorithm for solving integer and mixed integer optimization problems."
# Appl. Math. Comput. 212 (2009) 505-518

using Evolutionary, Test, Random

f1(x) = 2*x[1]+x[2]
lx = [0.0, 0]
ux = [1.6, 1]

cons(x) = [1.25-x[1]^2-x[2], x[1]+x[2]]
lc = [-Inf, -Inf]
uc = [0.0, 1.6]

tc = [Float64, Int]

cb = Evolutionary.ConstraintBounds(lx,ux,lc,uc)
c = MixedTypePenaltyConstraints(PenaltyConstraints(1e3, cb, cons), tc)
c = MixedTypePenaltyConstraints(WorstFitnessConstraints(cb, cons), tc)
init = ()->Real[rand(Float64), rand(0:1)]

opts = Evolutionary.Options(iterations=500, abstol=1e-5)
mthd = GA(populationSize=40, crossoverRate=0.8, mutationRate=0.05, selection=susinv, crossover=MILX(), mutation=MIPM(lx,ux))

Random.seed!(10);
result = Evolutionary.optimize(f1, c, init, mthd, opts)
@test minimum(result) ≈ 2.0 atol=1e-1
@test Evolutionary.minimizer(result) ≈ [0.5, 1] atol=1e-1

#-----------------------------------------------------------------

f8(z) = sum((z[4:end-1].-1).^2) - log(z[end]+1) + sum((z[i]-i)^2 for i in 1:3)

@test f8(Real[0.2, 1.280624, 1.954483, 1, 0, 0, 1]) ≈ 3.55746 atol=1e-6

cons(z) = [
    sum(z[1:end-1]),
    sum(v->v^2, z[[1,2,3,6]]),
    z[1]+z[4],
    z[2]+z[5],
    z[3]+z[6],
    z[1]+z[7],
    z[2]^2+z[5]^2,
    z[3]^2+z[6]^2,
    z[2]^2+z[6]^2,
]

lx = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
ux = [Inf, Inf, Inf, 1.0, 1.0, 1.0, 1.0]

lc = Float64[]
uc = [5.0, 5.5, 1.2, 1.8, 2.5, 1.2, 1.64, 4.25, 4.64]

tc = [Float64, Float64, Float64, Int, Int, Int, Int]

cb = Evolutionary.ConstraintBounds(lx,ux,lc,uc)
c = MixedTypePenaltyConstraints(PenaltyConstraints(1e2, cb, cons), tc)
c = MixedTypePenaltyConstraints(WorstFitnessConstraints(cb, cons), tc)
init = ()->Real[rand(Float64,3); rand(0:1,4)]

opts = Evolutionary.Options(iterations=2000)
mthd = GA(populationSize=100, ɛ=0.05, crossoverRate=0.8, mutationRate=0.01, selection=susinv, crossover=MILX(0.0,0.5,0.3), mutation=MIPM(lx,ux))

Random.seed!(235);
result = Evolutionary.optimize(f8, c, init, mthd, opts)
@test minimum(result) < 4

# [Evolutionary.optimize(f8, c, init, mthd, opts) |> minimum for i in 1:100] |> minimum
# uc .- (Evolutionary.minimizer(result) |> cons)
# Evolutionary.isfeasible(c.penalty, Evolutionary.minimizer(result))
# Evolutionary.minimizer(result) |> f8