text:
#
# variables.jl -
#
# Julia interface to Tcl variables.
#

# Manage to make any Tcl interpreter usable as an indexable collection with
# respect to its global variables.

Base.getindex(interp::TclInterp, name) = getvar(interp, name)

Base.getindex(interp::TclInterp, name1, name2) = getvar(interp, name1, name2)

Base.setindex!(interp::TclInterp, value, name) = begin
    setvar(interp, name, value)
    interp
end

Base.setindex!(interp::TclInterp, value, name1, name2) = begin
    setvar(interp, name1, name2, value)
    interp
end

Base.setindex!(interp::TclInterp, ::Nothing, name) = begin
    unsetvar(interp, name; nocomplain=true)
    interp
end

Base.setindex!(interp::TclInterp, ::Nothing, name1, name2) = begin
    unsetvar(interp, name1, name2; nocomplain=true)
    interp
end

Base.haskey(interp::TclInterp, name) = exists(interp, name)

Base.haskey(interp::TclInterp, name1, name2) = exists(interp, name1, name2)

const VARIABLE_FLAGS = TCL_GLOBAL_ONLY|TCL_LEAVE_ERR_MSG

"""
```julia
Tcl.getvar([T=Any,][interp,] name1 [,name2])
```

yields the value of the global variable `name1` or `name1(name2)` in Tcl
interpreter `interp` or in the initial interpreter if this argument is omitted.

Optional argument `T` (`Any` by default) can be used to specify the type of the
returned value.  Some possibilities are:

* If `T` is `Any`, the type of the returned value is determined so as to best
  reflect that of the Tcl variable.

* If `T` is `TclObj`, a managed Tcl object is returned.  This is the most
  efficient if the returned value is intended to be put in a Tcl list or to be
  an argument of a Tcl script or command.

* If `T` is `Bool`, a boolean value is returned.

* If `T` is `String`, a string is returned.

* If `T` is `Char`, a single character is returned (an exception is thrown if
  Tcl object is not a single character string).

* If `T <: Integer`, an integer value of type `T` is returned.

* If `T <: AbstractFloat`, a floating-point value of type `T` is returned.

* If `T` is `Vector`, a vector of values is returned (the Tcl object is
  converted into a list if necessary).

* `TclObj` to get a managed Tcl object;

* If `T` is `Vector`, a vector of values is returned (the Tcl object is
  converted into a list if necessary).

Note that, except if `T` is `Any` or `TclObj`, a conversion of the Tcl
object stored by the variable may be needed.

See also: [`Tcl.exists`](@ref), [`Tcl.setvar`](@ref), [`Tcl.unsetvar`](@ref).

"""
getvar(name::Name, args...) = getvar(getinterp(), name, args...)

getvar(::Type{T}, name::Name, args...) where {T} =
    getvar(T, getinterp(), name, args...)

getvar(interp::TclInterp, name::Name, args...) =
    getvar(Any, interp, name, args...)

function getvar(::Type{T}, interp::TclInterp, name::Name) where {T}
    ptr = __getvar(interp, name, VARIABLE_FLAGS)
    ptr != C_NULL || Tcl.error(interp)
    return __objptr_to(T, ptr)
end

function getvar(::Type{T}, interp::TclInterp,
                name1::Name, name2::Name) where {T}
    ptr = __getvar(interp, name1, name2, VARIABLE_FLAGS)
    ptr != C_NULL || Tcl.error(interp)
    return __objptr_to(T, ptr)
end

# Tcl_GetVar would yield an incorrect result if the variable value has embedded
# nulls and symbol Tcl_ObjGetVar2Ex does not exist in the library (despite what
# says the doc.).  Hence we always use Tcl_ObjGetVar2 which may require to
# temporarily convert the variable name into a string object.  There is no loss
# of performances as it turns out that Tcl_GetVar, Tcl_GetVar2 and
# Tcl_ObjGetVar2Ex all call Tcl_ObjGetVar2.
#
# Tcl_ObjGetVar2 does not manage the reference count of the variable name
# parts.
#
# Since all arguments are passed as pointers to Tcl object, we have to take
# care of correctly unreference temporary objects.  Since there can be no
# failure for argument types allowed for variable name parts, we can avoid the
# overhead of the `try ... catch ... finally` statements.

function __getvar(interp::TclInterp, name::Name, flags::Integer)
    nameptr = Tcl_IncrRefCount(__objptr(name))
    objptr = Tcl_ObjGetVar2(interp.ptr, nameptr, C_NULL, flags)
    Tcl_DecrRefCount(nameptr)
    return objptr
end

function __getvar(interp::TclInterp, name1::Name, name2::Name, flags::Integer)
    name1ptr = Tcl_IncrRefCount(__objptr(name1))
    name2ptr = Tcl_IncrRefCount(__objptr(name2))
    objptr = Tcl_ObjGetVar2(interp.ptr, name1ptr, name2ptr, flags)
    Tcl_DecrRefCount(name1ptr)
    Tcl_DecrRefCount(name2ptr)
    return objptr
end

"""
```julia
Tcl.setvar([interp,] name1, [name2,] value)
```

set global variable `name1` or `name1(name2)` to be `value` in Tcl interpreter
`interp` or in the initial interpreter if this argument is omitted.  The result
is `nothing`.

See [`Tcl.getvar`](@ref) for details about allowed variable names.

See also: [`Tcl.getvar`](@ref), [`Tcl.exists`](@ref), [`Tcl.unsetvar`](@ref).

"""
setvar(name::Name, args...) = setvar(getinterp(), name, args...)

function setvar(interp::TclInterp, name::Name, value)
    if __setvar(interp, name, value, VARIABLE_FLAGS) == C_NULL
        Tcl.error(interp)
    end
    return nothing
end

function setvar(interp::TclInterp, name1::Name, name2::Name, value)
    if __setvar(interp, name1, name2, value, VARIABLE_FLAGS) == C_NULL
        Tcl.error(interp)
    end
    return nothing
end

# Like Tcl_ObjGetVar2Ex, Tcl_ObjSetVar2Ex may be not found in library so we
# avoid using it.  In fact, it turns out that Tcl_SetVar, Tcl_SetVar2 and
# Tcl_ObjSetVar2Ex call Tcl_ObjGetVar2 to do their stuff, so we only use
# Tcl_ObjGetVar2 with no loss of performances as it turns out that Tcl_SetVar
# calls Tcl_ObjGetVar2.
#
# Tcl_ObjSetVar2 does not manage the reference count of the variable name
# parts but do manage the reference count of the variable value.
#
# Same remarks as for `__getvar` about correctly unreferencing temporary
# objects and avoiding `try ... catch` statements.

function __setvar(interp::TclInterp, name::Name, value, flags::Integer)
    valueptr = Tcl_IncrRefCount(__objptr(value)) # must be first
    nameptr = Tcl_IncrRefCount(__objptr(name))
    objptr = Tcl_ObjSetVar2(interp.ptr, nameptr, C_NULL, valueptr, flags)
    Tcl_DecrRefCount(valueptr)
    Tcl_DecrRefCount(nameptr)
    return objptr
end

function __setvar(interp::TclInterp, name1::Name, name2::Name,
                  value, flags::Integer)
    valueptr = Tcl_IncrRefCount(__objptr(value)) # must be first
    name1ptr = Tcl_IncrRefCount(__objptr(name1))
    name2ptr = Tcl_IncrRefCount(__objptr(name2))
    objptr = Tcl_ObjSetVar2(interp.ptr, name1ptr, name2ptr, valueptr, flags)
    Tcl_DecrRefCount(valueptr)
    Tcl_DecrRefCount(name1ptr)
    Tcl_DecrRefCount(name2ptr)
    return objptr
end

"""
```julia
Tcl.unsetvar([interp,] name1 [,name2]; nocomplain=false)
```

deletes global variable `name1` or `name1(name2)` in Tcl interpreter `interp`
or in the initial interpreter if this argument is omitted.

Keyword `nocomplain` can be set true to ignore errors.

See also: [`Tcl.getvar`](@ref), [`Tcl.exists`](@ref), [`Tcl.setvar`](@ref).

"""
unsetvar(args...; kwds...) = unsetvar(getinterp(), args...; kwds...)

function unsetvar(interp::TclInterp, name::Name; nocomplain::Bool=false)
    flags = (nocomplain ? TCL_GLOBAL_ONLY : (TCL_GLOBAL_ONLY|TCL_LEAVE_ERR_MSG))
    status = __unsetvar(interp, __string(name), flags)
    if status != TCL_OK && ! nocomplain
        Tcl.error(interp)
    end
    return nothing
end

function unsetvar(interp::TclInterp, name1::Name, name2::Name;
                  nocomplain::Bool=false)
    flags = (nocomplain ? TCL_GLOBAL_ONLY : (TCL_GLOBAL_ONLY|TCL_LEAVE_ERR_MSG))
    status = __unsetvar(interp, __string(name1), __string(name2), flags)
    if status != TCL_OK && ! nocomplain
        Tcl.error(interp)
    end
    return nothing
end

# `TclUnsetVarObj2` would be the function to call here but, unfortunately, only
# `Tcl_UnsetVar` and `Tcl_UnsetVar2` are available which both require strings
# for the variable name parts.

function __unsetvar(interp::TclInterp, name::String, flags::Integer)
    if (ptr = __cstring(name)[1]) != C_NULL
        status = Tcl_UnsetVar(interp.ptr, ptr, flags)
    else
        status = Tcl.eval(TclStatus, interp, "unset {$name}")
    end
    return status
end

function __unsetvar(interp::TclInterp, name1::String, name2::String,
                    flags::Integer)
    if ((ptr1 = __cstring(name1)[1]) != C_NULL &&
        (ptr2 = __cstring(name2)[1]) != C_NULL)
        status = Tcl_UnsetVar2(interp.ptr, ptr1, ptr2, flags)
    else
        status = Tcl.eval(TclStatus, "unset {$name1($name2)}")
    end
    return status
end

"""
`julia
Tcl.exists([interp,] name1 [, name2])
```

checks whether global variable `name1` or `name1(name2)` is defined in Tcl
interpreter `interp` or in the initial interpreter if this argument is omitted.

See also: [`Tcl.getvar`](@ref), [`Tcl.setvar`](@ref), [`Tcl.unsetvar`](@ref).

"""
exists(args...) = exists(getinterp(), args...)

function exists(interp::TclInterp, name::Name)
    return (__getvar(interp, name, TCL_GLOBAL_ONLY) != C_NULL)
end

function exists(interp::TclInterp, name1::Name, name2::Name)
    return (__getvar(interp, name1, name2, TCL_GLOBAL_ONLY) != C_NULL)
end

"""
```julia
__string(str)
```

yields a `String` instance of `str`.

"""
__string(str::String) = str
__string(str::Union{AbstractString,Symbol,TclObj{String}}) = string(str)

"""
```julia
__cstring(str) -> ptr, siz
```

checks whether `str` is a valid C-string (i.e., has no embedded nulls)
and yields its base address and size.  If `str` is not eligible as a C string,
`(Ptr{Cchar}(0), 0)` is returned.

"""
function __cstring(str::String)
    ptr, siz = Base.unsafe_convert(Ptr{Cchar}, str), sizeof(str)
    if Base.containsnul(ptr, siz)
        ptr, siz = Ptr{Cchar}(0), zero(siz)
    end
    return ptr, siz
end

__cstring(sym::Symbol) = __cstring(string(sym))


text:
module CurrenciesBase

using FixedPointDecimals

import Base: +, -, *, /, ==

# Exports
export AbstractMonetary, Monetary
export currency, decimals, majorunit, @usingcurrencies
export currencyinfo, iso4217num, iso4217alpha, shortsymbol, longsymbol
export newcurrency!, @usingcustomcurrency

# Currency data
include("data/currencies.jl")
include("data/locale.jl")
include("data/symbols.jl")

# Monetary type, currencies, and arithmetic
include("monetary.jl")
include("currency.jl")
include("arithmetic.jl")
include("mixed.jl")

# Custom currencies and macros
include("usingcurrencies.jl")
include("custom.jl")

# Display
include("formatting.jl")

end  # module CurrenciesBase

ext: jl


text:
struct MS_PSO{T,S,fType} <: Optimizer

    # Optimization problem
    prob::Problem{fType,S}

    # Vector of particle swarms
    swarmVec::Vector{Swarm{T}}

    # Settings
    initMethod::Symbol
    updateMethod::Symbol

    # PSO specific parameters
    inertiaRange::Tuple{T,T}
    minNeighborFrac::T
    selfAdjustWeight::T
    socialAdjustWeight::T

    # Results data
    results::Results{T}

    function MS_PSO{T}(prob::Problem{fType,S}, numParticlesPerSwarm::Integer,
        numSwarms::Integer, initMethod::Symbol, updateMethod::Symbol,
        inertiaRange::Tuple{T,T}, minNeighborFrac::T, selfAdjustWeight::T,
        socialAdjustWeight::T) where {T,S,fType}

        # Initialize Swarm vector
        N = length(prob.LB)
        swarmVec = [Swarm{T}(N, numParticlesPerSwarm) for i in 1:numSwarms]

        # Initialize results
        results = Results{T}(undef, N)

        # Check that init and update methods exist
        if initMethod != :Uniform && initMethod != :LogisticsMap
            throw(ArgumentError("$initMethod is not a PSO initialization method."))
        end
        if updateMethod != :MATLAB
            throw(ArgumentError("$updateMethod is not a PSO update method."))
        end

        return new{T,S,fType}(prob, swarmVec, initMethod, updateMethod,
                        inertiaRange, minNeighborFrac, selfAdjustWeight,
                        socialAdjustWeight, results)
    end
end

# ===== Constructors

function MS_PSO(prob::Problem{fType,S}; numParticlesPerSwarm = 100, numSwarms = 4,
    inertiaRange = (0.1, 1.1), minNeighborFrac = 0.25, selfAdjustWeight = 1.49,
    socialAdjustWeight = 1.49, initMethod::Symbol = :Uniform,
    updateMethod::Symbol = :MATLAB) where {S,fType}

    # Error checking
    length(inertiaRange) == 2 || throw(ArgumentError("inertiaRange must be of length 2."))
    minNeighborFrac > 0       || throw(ArgumentError("minNeighborFrac must be > 0."))

    # Type info
    T = typeof(1.0)
    nIRange = (T(inertiaRange[1]), T(inertiaRange[2]))

    # Call constructor
    return MS_PSO{T}(prob, numParticlesPerSwarm, numSwarms, initMethod, updateMethod,
        nIRange, T(minNeighborFrac), T(selfAdjustWeight), T(socialAdjustWeight))
end

# ===== Methods

function _optimize!(mspso::MS_PSO, opts::Options)

    # Initialize PSO
    initialize!(mspso, opts)

    # Iterate PSO
    iterate!(mspso, opts)

    # Display results and return
    if opts.display
        display(mspso.results)
    end

    return mspso.results
end

function initialize!(mspso::MS_PSO, opts::Options)
    @inbounds begin
    # Initialize position and velocities
    if mspso.initMethod == :Uniform
        for i in 1:length(mspso.swarmVec)
            uniformInitialization!(mspso.swarmVec[i], mspso.prob, opts)
        end
    elseif mspso.initMethod == :LogisticsMap
        for i in 1:length(mspso.swarmVec)
            logisticsMapInitialization!(mspso.swarmVec[i], mspso.prob, opts)
        end
    else
        throw(ArgumentError("PSO initialization method is not implemented."))
    end

    # Iterate through each swarm
    # Note: This implementation is not entirely optimal as we will
    # be waiting for each swarm to finish initializing before starting the next
    for i in 1:length(mspso.swarmVec)
        # Evaluate Objective Function
        feval!(mspso.swarmVec[i], mspso.prob.f, opts; init = true)

        # Set swarm global best obj. func. value to Inf
        mspso.swarmVec[i].b = Inf

        # Set global best
        setGlobalBest!(mspso.swarmVec[i])

        # Initialize neighborhood size
        mspso.swarmVec[i].n = max(2, floor(length(mspso.swarmVec[i])*mspso.minNeighborFrac))

        # Initial inertia
        if mspso.inertiaRange[2] > 0
            mspso.swarmVec[i].w = (mspso.inertiaRange[2] > mspso.inertiaRange[1]) ?
                                mspso.inertiaRange[2] : mspso.inertiaRange[1]
        else
            mspso.swarmVec[i].w = (mspso.inertiaRange[2] < mspso.inertiaRange[1]) ?
                                mspso.inertiaRange[2] : mspso.inertiaRange[1]
        end

        # Initialize self and social adjustment
        mspso.swarmVec[i].y₁ = mspso.selfAdjustWeight
        mspso.swarmVec[i].y₂ = mspso.socialAdjustWeight
    end

    # Print Status # NEED TO IMPLEMENT A NEW printStatus for Vector{Swarm}
    if opts.display
        printStatus(mspso.swarmVec, 0.0, 0, 0)
    end

    end
    return nothing
end

function iterate!(mspso::MS_PSO, opts::Options)

    # Initialize time and iteration counter
    t0  = time()
    stallT0 = t0*ones(length(mspso.swarmVec))
    iters = 0
    stallIters = fill!(Vector{Int}(undef, length(mspso.swarmVec)), 0)
    hasStalled = BitVector(undef, length(mspso.swarmVec))
    fStall = Inf*ones(length(mspso.swarmVec))

    # Compute minimum neighborhood size
    minNeighborSize = max(2, floor(length(mspso.swarmVec[1])*mspso.minNeighborFrac))

    # Begin loop
    exitFlag = 0
    while exitFlag == 0

        # Update iteration counter
        iters += 1

        # Prepare all swarms for evaluation
        for i in 1:length(mspso.swarmVec)
            # Update swarm velocities
            updateVelocities!(mspso.swarmVec[i])

            # Update positions
            step!(mspso.swarmVec[i])

            # Enforce Bounds
            if !(all(isinf.(mspso.prob.LB)) && all(isinf.(mspso.prob.UB)))
                enforceBounds!(mspso.swarmVec[i], mspso.prob.LB, mspso.prob.UB)
            end
        end

        # Evaluate objective function for all swarms
        # Will need to experiment with multithreadding here!
        for i in 1:length(mspso.swarmVec)
            feval!(mspso.swarmVec[i], mspso.prob.f, opts)
        end

        # Update all swarms
        for i in 1:length(mspso.swarmVec)
            # Update global best
            flag::Bool = setGlobalBest!(mspso.swarmVec[i])

            # Update stall counter and neighborhood
            if flag
                mspso.swarmVec[i].c = max(0, mspso.swarmVec[i].c - 1)
                mspso.swarmVec[i].n = minNeighborSize
            else
                mspso.swarmVec[i].c += 1
                mspso.swarmVec[i].n = min(mspso.swarmVec[i].n + minNeighborSize,
                                      length(mspso.swarmVec[i]) - 1)
            end

            # Update inertia
            if mspso.swarmVec[i].c < 2
                mspso.swarmVec[i].w *= 2.0
            elseif mspso.swarmVec[i].c > 5
                mspso.swarmVec[i].w /= 2.0
            end

            # Ensure inertia is in inertia range
            if  mspso.swarmVec[i].w > mspso.inertiaRange[2]
                mspso.swarmVec[i].w = mspso.inertiaRange[2]
            elseif mspso.swarmVec[i].w < mspso.inertiaRange[1]
                mspso.swarmVec[i].w = mspso.inertiaRange[1]
            end

            # Track stalling
            if fStall[i] - mspso.swarmVec[i].b > opts.funcTol
                fStall[i] = mspso.swarmVec[i].b
                stallIters[i] = 0
                stallT0[i] = time()
            else
                stallIters[i] += 1
            end

            # Check for stalling
            if stallIters[i] >= opts.maxStallIters
                hasStalled[i] = true
                interMingle!(mspso.swarmVec, i)
            end

       end

        # Stopping criteria
        if all(hasStalled)
            exitFlag = 1
        elseif iters >= opts.maxIters
            exitFlag = 2
        elseif any([mspso.swarmVec[i].b <= opts.objLimit for i in 1:length(mspso.swarmVec)])
            exitFlag = 3
        elseif all([time() - stallT0[i] >= opts.maxStallTime for i in 1:length(mspso.swarmVec)])
            exitFlag = 4
        elseif time() - t0 >= opts.maxTime
            exitFlag = 5
        end

        # Output Status
        if opts.display && iters % opts.displayInterval == 0
            printStatus(mspso.swarmVec, time() - t0, iters, stallIters)
        end
    end

    # Set results
    setResults!(mspso, iters, time() - t0, exitFlag)

    return nothing
end

# Could intermingle stalled swarm with best swarm or random swarm. Choosing random swarm for now!
function interMingle!(swarmVec::Vector{Swarm{T}}, stallIdx) where {T}
    # Choosing swarm to intermingle with stalled swarm
    swarmChosen = false
    swarmIdx = stallIdx
    while !swarmChosen
        swarmIdx = rand(1:length(swarmVec))
        if swarmIdx != stallIdx
            swarmChosen = true
        end
    end

    # Intermingle small subset of chosen swarm (10%)
    n = length(swarmVec[1])
    nSubSet = floor(0.1*n)
    for i in 1:nSubSet
        idx = rand(1:n)
        swarmVec[stallIdx][idx].p .= swarmVec[swarmIdx][idx].p
        swarmVec[stallIdx][idx].fp = swarmVec[swarmIdx][idx].fp
    end

    return nothing
end

function printStatus(swarmVec::Vector{Swarm{T}}, time, iter, stallCount::Vector{Int}) where {T}
    bestF = minimum([swarmVec[i].b for i in 1:length(swarmVec)])
    fspec1 = FormatExpr("Time Elapsed: {1:f} sec, Iteration Number: {2:d}, Function Evaluations: {3:d}")
    fspec2 = FormatExpr("Max Stall Iterations: {1:d}, Global Best: {2:e}")
    printfmtln(fspec1, time, iter, (iter + 1)*length(swarmVec[1])*length(swarmVec))
    printfmtln(fspec2, maximum([stallCount[i] for i in 1:length(swarmVec)]), bestF)
    println(" ")
end

function printStatus(swarmVec::Vector{Swarm{T}}, time, iter, stallCount::Int) where {T}
    bestF = minimum([swarmVec[i].b for i in 1:length(swarmVec)])
    fspec1 = FormatExpr("Time Elapsed: {1:f} sec, Iteration Number: {2:d}, Function Evaluations: {3:d}")
    fspec2 = FormatExpr("Max Stall Iterations: {1:d}, Global Best: {2:e}")
    printfmtln(fspec1, time, iter, (iter + 1)*length(swarmVec[1])*length(swarmVec))
    printfmtln(fspec2, stallCount, bestF)
    println(" ")
end

function setResults!(mspso::MS_PSO, iters::Int, Δt::AbstractFloat, exitFlag::Int)
    bestIdx = -1
    bestF = Inf
    for i in 1:length(mspso.swarmVec)
        if mspso.swarmVec[i].b < bestF
            bestF = mspso.swarmVec[i].b
            bestIdx = i
        end
    end
    mspso.results.fbest = bestF
    mspso.results.xbest .= mspso.swarmVec[bestIdx].d
    mspso.results.iters = iters
    mspso.results.time = Δt
    mspso.results.exitFlag = exitFlag

    return nothing
end


text:
function is_integer_feasible(m::GlobalOptimizer)
    bool = true
    atol = _integer_abs_tol(m)
    rtol = _integer_rel_tol(m)
    for i = 1:_variable_num(BranchVar(), m)
        if is_integer(BranchVar(), m, i)
            xsol = _lower_solution(BranchVar(), m, i)
            if isapprox(floor(xsol), xsol; atol = atol, rtol = rtol)
                continue
            elseif isapprox(ceil(xsol), xsol; atol = atol, rtol = rtol)
                continue
            else
                bool &= false
                break
            end
        end
    end
    return bool
end

"""
$(SIGNATURES)

Default check to see if the upper bounding problem should be run. By default,
The upper bounding problem is run on every node up to depth `upper_bounding_depth`
and is triggered with a probability of `0.5^(depth - upper_bounding_depth)`
afterwards.
"""
function default_nlp_heurestic(m::GlobalOptimizer)
    bool = false
    ubd_limit = m._parameters.upper_bounding_depth
    n = _current_node(m)
    if is_integer_feasible(m)
        Δdepth = n.depth - n.cont_depth
        bool |= (Δdepth <= ubd_limit)
        bool |= (rand() < 0.5^(Δdepth - ubd_limit))
    end
    return bool
end

"""
$(SIGNATURES)

Default upper bounding problem which simply calls `solve_local_nlp!` to solve
the nlp locally.
"""
function upper_problem!(t::ExtensionType, m::GlobalOptimizer)
    if !default_nlp_heurestic(m)
        m._upper_feasibility = false
        m._upper_objective_value = Inf
    else
        solve_local_nlp!(m)
    end
    return
end
upper_problem!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = upper_problem!(_ext_typ(m), m)

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Julia script for eEOM demo
# Author: Yen-Lin Chen

# Load the dependencies
using Plots, DelimitedFiles, StatsBase, JLD
include("ga.jl");

# Load data
pool1, pool2 = load("pool.jld", "pool1", "pool2");
data1, data2 = readdlm("data1.dat"), readdlm("data2.dat");
q = collect(0.0:0.005:1.25);

## First look at the [data1, pool1] set

# Plot the data and see what's going on
p = plot(size=(400, 300), dpi=600, grid=false, legend=:outerbottomright, framestyle=:box);
s = 0.4;
for i = 1:14
    plot!(data1[:,1], s^(i-1) * data1[:,2], yscale=:log10, size=(800,900), lw=2.0, c=:black, linealpha=0.5, lab="");
    plot!(q, s^(i-1)*pool1[:,7i-6], lw=1.5, yscale=:log10, lab="Conf. $(12+i)");
end
xlims!(0.05, 0.8)
ylims!(0.8e-6, 10)
ylabel!("I(q) (A.U.)")
xlabel!("q (1/A)")

# Run the eEOM structural refinement
idx = findall(data1[:,1] .<= 0.8);
tmp = data1[idx, :];
ensize = 10;
cycles = 50;
survivors1 = zeros(Int64, ensize, cycles);
fitness1 = zeros(cycles);

t = @elapsed for i = 1:cycles
    @info("===== Cycle number $i ====");
    fitness1[i], survivor = Refine(tmp, [q pool1], ensize, 0.15; population=10, maxiter=100000, reports=10);
    survivors1[:,i] = survivor[:,1];
end

println("- eEOM for $cycles cycles and $(size(pool1, 2)) conformaitons took $(round(t, digits=2)) seconds.");

# Plot the eEOM results
p = scatter(data1[:,1], data1[:,2], yerror=data1[:,3], yscale=:log10, markersize=5.0, c=:black, markeralpha=0.6, linealpha=0.25, lab="Experiment", size=(300,450), legendfont=font(16), titlefont=font(20));
plot!(q, mean(pool1[:, survivors1[:, findall(fitness1 .== minimum(fitness1))[1]]], dims=2), lw=5, linealpha=0.8, lab="eEOM Fit")
plot!(size=(600, 800), dpi=600, grid=false, legend=true, framestyle=:box);
xlims!(0.05, 0.8);
ylims!(0.1, 10)
xlabel!("q (1/A)");
ylabel!("I(q) (A.U.)")

# Plot the best ensemble
best = survivors1[:, findall(fitness1 .== minimum(fitness1))[1]];
dic = countmap(reshape(best, (length(best), )));
b = bar(collect(keys(dic)), collect(values(dic)) / length(best), title="Cluster Selection: Best", xticks=collect(1:2:98), bar_width=0.5);
plot!(size=(500, 300), dpi=600, grid=false, legend=false, framestyle=:box);
ylabel!("Population Fraction")

# Plot the all ensemble
dic = countmap(reshape(survivors1, (length(survivors1), )));
b = bar(collect(keys(dic)), collect(values(dic)) / length(best), title="Cluster Selection: All", xticks=collect(1:2:98), bar_width=0.5);
plot!(size=(500, 300), dpi=600, grid=false, legend=false, framestyle=:box);
ylabel!("Population Fraction")

# And then do whatever you like from the results

## Second look at the [data2, pool2] set

# Plot the data and see what's going on
p = plot(size=(400, 300), dpi=600, grid=false, legend=:outerbottomright, framestyle=:box);
s = 0.4;
for i = 1:14
    plot!(data2[:,1], s^(i-1) * data2[:,2], yscale=:log10, size=(800,900), lw=2.0, c=:black, linealpha=0.5, lab="");
    plot!(q, s^(i-1)*pool2[:,7i-6], lw=1.5, yscale=:log10, lab="Conf. $(12+i)");
end
xlims!(0.05, 0.8)
ylims!(0.8e-6, 10)
ylabel!("I(q) (A.U.)")
xlabel!("q (1/A)")

# Run the eEOM structural refinement
idx = findall(data2[:,1] .<= 1.0);
tmp = data2[idx, :];
ensize = 10;
cycles = 50;
survivors2 = zeros(Int64, ensize, cycles);
fitness2 = zeros(cycles);

t = @elapsed for i = 1:cycles
    @info("===== Cycle number $i ====");
    fitness2[i], survivor = Refine(tmp, [q pool2], ensize, 0.10; population=10, maxiter=100000, reports=10);
    survivors2[:,i] = survivor[:,1];
end

println("- eEOM for $cycles cycles and $(size(pool2, 2)) conformaitons took $(round(t, digits=2)) seconds.");

# Plot the eEOM results
p = scatter(data2[:,1], data2[:,2], yerror=data2[:,3], yscale=:log10, markersize=5.0, c=:black, markeralpha=0.6, linealpha=0.25, lab="Experiment", size=(300,450), legendfont=font(16), titlefont=font(20));
plot!(q, mean(pool2[:, survivors2[:, findall(fitness2 .== minimum(fitness2))[1]]], dims=2), lw=5, linealpha=0.8, lab="eEOM Fit")
plot!(size=(600, 800), dpi=600, grid=false, legend=true, framestyle=:box);
xlims!(0.05, 0.8);
ylims!(0.1, 10)
xlabel!("q (1/A)");
ylabel!("I(q) (A.U.)")

# Plot the best ensemble
best = survivors2[:, findall(fitness2 .== minimum(fitness2))[1]];
dic = countmap(reshape(best, (length(best), )));
b = bar(collect(keys(dic)), collect(values(dic)) / length(best), title="Cluster Selection: Best", xticks=collect(1:4:98), bar_width=0.5);
plot!(size=(500, 300), dpi=600, grid=false, legend=false, framestyle=:box);
ylabel!("Population Fraction")

# Plot the all ensemble
dic = countmap(reshape(survivors2, (length(survivors2), )));
b = bar(collect(keys(dic)), collect(values(dic)) / length(best), title="Cluster Selection: All", xticks=collect(1:4:98), bar_width=0.5);
plot!(size=(500, 300), dpi=600, grid=false, legend=false, framestyle=:box);
ylabel!("Population Fraction")

# And then do whatever you like from the results

# EOF


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    getsourcenode, getsourceedge, getunionnode, getunionedge,
    disjointunion, disjointunion!

struct DisjointUnionNode <: AbstractNode
    source::Int
    sourcekey::Int
end

struct DisjointUnionEdge <: UndirectedEdge
    source::Int
    sourcekey::Int
end

struct DisjointUnionGraph <: OrderedGraph
    neighbormap::Vector{Dict{Int,Int}}
    edges::Vector{Tuple{Int,Int}}
    nodeattrs::Vector{DisjointUnionNode}
    edgeattrs::Vector{DisjointUnionEdge}
    cache::Dict{Symbol,Any}
end

getsourcenode(graph, i) = graph.nodeattrs[i]
getsourceedge(graph, i) = graph.edgeattrs[i]

function getunionnode(graph::DisjointUnionGraph, source, key)
    for (i, node) in enumerate(nodeattrs(graph))
        if node.source == source && node.sourcekey == key
            return i
        end
    end
end

function getunionedge(graph::DisjointUnionGraph, source, key)
    for (i, edge) in enumerate(edgeattrs(graph))
        if edge.source == source && edge.sourcekey == key
            return i
        end
    end
end

"""
    disjointunion(g1::UndirectedGraph, g2::UndirectedGraph,
        G::UndirectedGraph...) -> DisjointUnionGraph

Generate disjoint union graph of given graphs. The new graph with type
`DisjointUnionGraph` retains mapping to the original graphs as nodes and edges
attributes.
"""
function disjointunion(
            g1::UndirectedGraph, g2::UndirectedGraph, G::UndirectedGraph...)
    U = DisjointUnionGraph([], [], [], [], Dict())
    for (src, g) in enumerate([g1, g2, G...])
        nsort = sort(collect(nodeset(g)))
        esort = sort(collect(edgeset(g)))
        nmap = Dict(n => i + nodecount(U) for (i, n) in enumerate(nsort))
        emap = Dict(e => i + edgecount(U) for (i, e) in enumerate(esort))
        nbrs = [neighbors(g, n) for n in nsort]
        edges = [getedge(g, e) for e in esort]
        append!(U.neighbormap, relabel(nbrs, nmap, emap))
        append!(U.edges, relabel(edges, nmap))
        append!(U.nodeattrs, [DisjointUnionNode(src, k) for k in nsort])
        append!(U.edgeattrs, [DisjointUnionEdge(src, k) for k in esort])
    end
    return U
end

"""
    disjointunion!(g1::T, g2::T, G::T...) where {T<:OrderedGraph} -> T

Generate disjoint union graph of given graphs. `g1` will be overwritten by the
union graph. Unlike non-destructive `disjointunion`, `g1` does not retain any
information about other given graphs but a bit faster.
"""
function disjointunion!(g1::T, g2::T, G::T...) where {T<:OrderedGraph}
    for g in [g2, G...]
        noff = nodecount(g1)
        nmap = Dict(i => i + noff for i in 1:nodecount(g))
        eoff = edgecount(g1)
        emap = Dict(i => i + eoff for i in 1:edgecount(g))
        append!(g1.neighbormap, relabel(g.neighbormap, nmap, emap))
        append!(g1.edges, relabel(g.edges, nmap))
        isdefined(g1, :nodeattrs) && append!(g1.nodeattrs, g.nodeattrs)
        isdefined(g1, :edgeattrs) && append!(g1.edgeattrs, g.edgeattrs)
    end
end

function relabel(edges::Vector{Tuple{Int,Int}}, nmap::Dict{Int,Int})
    es = Tuple{Int,Int}[]
    for (u, v) in edges
        push!(es, (nmap[u], nmap[v]))
    end
    return es
end

function relabel(neighbormap::Vector{Dict{Int,Int}},
                 nmap::Dict{Int,Int}, emap::Dict{Int,Int})
    ns = Dict{Int,Int}[]
    for nbr in neighbormap
        d = Dict{Int,Int}()
        for (inc, adj) in nbr
            d[emap[inc]] = nmap[adj]
        end
        push!(ns, d)
    end
    return ns
end


text:
@testset "misc Permutation functionality" begin
    G = PermGroup([perm"(1,2)", perm"(1,2,3,4)"])
    S = gens(G)

    @test permtype(S[1]) == [2,1,1]
    @test permtype(S[2]) == [4]
    @test sign(S[1]) == -1
    @test sign(S[2]) == -1
    @test sign(S[2]^2) == 1
end


text:
"""
    dataset(name::String)::DataFrame

Load the dataset associated with `name`.

Some datasets available:
 - `50408`: Gauging (water levels and corresponding discharges) of the Sainte-Anne river in the province of Quebec (Canada).

# Examples
```julia-repl
julia> RatingCurves.dataset("50408")
```
"""
function dataset(name::String)::DataFrame

    filename = joinpath(dirname(@__FILE__), "..", "data", string(name, ".csv"))
    if isfile(filename)
        # return DataFrame!(CSV.File(filename))
        return CSV.read(filename, DataFrame)
    end
    error("There is no dataset with the name '$name'")

end


text:
"Explicit partitioned Runge-Kutta integrator."
immutable IntegratorEPRK{T} <: Integrator{T}
    equation::PODE{T}
    tableau::TableauEPRK{T}
    Δt::T

    q::Array{T,1}
    p::Array{T,1}
    y::Array{T,1}
    z::Array{T,1}
    Q::Array{T,2}
    P::Array{T,2}
    Y::Array{T,2}
    Z::Array{T,2}
    F::Array{T,2}
    G::Array{T,2}
    tQ::Array{T,1}
    tP::Array{T,1}
    tF::Array{T,1}
    tG::Array{T,1}

    function IntegratorEPRK(equation, tableau, Δt)
        D = equation.d
        S = tableau.s
        new(equation, tableau, Δt,
            zeros(T,D), zeros(T,D),
            zeros(T,D), zeros(T,D),
            zeros(T,D,S), zeros(T,D,S),
            zeros(T,D,S), zeros(T,D,S),
            zeros(T,D,S), zeros(T,D,S),
            zeros(T,D), zeros(T,D),
            zeros(T,D), zeros(T,D))
    end
end

function IntegratorEPRK{T}(equation::PODE{T}, tableau::TableauEPRK{T}, Δt::T)
    IntegratorEPRK{T}(equation, tableau, Δt)
end

"Compute Q stages of explicit partitioned Runge-Kutta methods."
function computeStageQ!(int::IntegratorEPRK, i::Int, jmax::Int, t)
    for j in 1:jmax
        for k in 1:int.equation.d
            int.Y[k,i] += int.tableau.a_q[i,j] * int.F[k,j]
        end
    end
    for k in 1:int.equation.d
        int.Q[k,i] = int.q[k] + int.Δt * int.Y[k,i]
    end
    simd_copy_xy_first!(int.tQ, int.Q, i)
    int.equation.g(t, int.tQ, int.tG)
    simd_copy_yx_first!(int.tG, int.G, i)
    nothing
end

"Compute P stages of explicit partitioned Runge-Kutta methods."
function computeStageP!(int::IntegratorEPRK, i::Int, jmax::Int, t)
    for j in 1:jmax
        for k in 1:int.equation.d
            int.Z[k,i] += int.tableau.a_p[i,j] * int.G[k,j]
        end
    end
    for k in 1:int.equation.d
        int.P[k,i] = int.p[k] + int.Δt * int.Z[k,i]
    end
    simd_copy_xy_first!(int.tP, int.P, i)
    int.equation.f(t, int.tP, int.tF)
    simd_copy_yx_first!(int.tF, int.F, i)
    nothing
end

"Integrate partitioned ODE with explicit partitioned Runge-Kutta integrator."
function integrate!{T}(int::IntegratorEPRK{T}, sol::SolutionPODE{T})
    # loop over initial conditions
    for m in 1:sol.n0
        local j::Int
        local tqᵢ::T
        local tpᵢ::T

        # copy initial conditions from solution
        for i in 1:sol.nd
            int.q[i] = sol[1, i, 0, m]
            int.p[i] = sol[2, i, 0, m]
        end

        for n in 1:sol.ntime
            # compute internal stages
            fill!(int.Y, zero(T))
            fill!(int.Z, zero(T))
            for i in 1:int.tableau.s
                tqᵢ = sol.t[n] + int.Δt * int.tableau.c_q[i]
                tpᵢ = sol.t[n] + int.Δt * int.tableau.c_p[i]

                if int.tableau.a_q[i,i] ≠ 0. && int.tableau.a_p[i,i] ≠ 0.
                    error("This is an implicit method!")
                elseif int.tableau.a_q[i,i] ≠ 0.
                    computeStageP!(int, i, i-1, tpᵢ)
                    computeStageQ!(int, i, i, tqᵢ)
                elseif int.tableau.a_p[i,i] ≠ 0.
                    computeStageQ!(int, i, i-1, tqᵢ)
                    computeStageP!(int, i, i, tpᵢ)
                else
                    computeStageQ!(int, i, i-1, tqᵢ)
                    computeStageP!(int, i, i-1, tpᵢ)
                end
            end

            # compute final update
            simd_mult!(int.y, int.F, int.tableau.b_q)
            simd_axpy!(int.Δt, int.y, int.q)

            simd_mult!(int.z, int.G, int.tableau.b_p)
            simd_axpy!(int.Δt, int.z, int.p)

            # copy to solution
            if mod(n, sol.nsave) == 0
                j = div(n, sol.nsave)
                for i in 1:sol.nd
                    sol[1, i, j, m] = int.q[i]
                    sol[2, i, j, m] = int.p[i]
                end
            end
        end
    end
    nothing
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Transfer <: Service
        name::String
        available::Bool
        requirement::Float64
        ext::Dict{String, Any}
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

# Arguments
- `name::String`
- `available::Bool`
- `requirement::Float64`
- `ext::Dict{String, Any}`
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct Transfer <: Service
    name::String
    available::Bool
    requirement::Float64
    ext::Dict{String, Any}
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function Transfer(name, available, requirement, ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    Transfer(name, available, requirement, ext, forecasts, InfrastructureSystemsInternal(), )
end

function Transfer(; name, available, requirement, ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    Transfer(name, available, requirement, ext, forecasts, )
end

# Constructor for demo purposes; non-functional.
function Transfer(::Nothing)
    Transfer(;
        name="init",
        available=false,
        requirement=0.0,
        ext=Dict{String, Any}(),
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

InfrastructureSystems.get_name(value::Transfer) = value.name
"""Get Transfer available."""
get_available(value::Transfer) = value.available
"""Get Transfer requirement."""
get_requirement(value::Transfer) = value.requirement
"""Get Transfer ext."""
get_ext(value::Transfer) = value.ext

InfrastructureSystems.get_forecasts(value::Transfer) = value.forecasts
"""Get Transfer internal."""
get_internal(value::Transfer) = value.internal

InfrastructureSystems.set_name!(value::Transfer, val::String) = value.name = val
"""Set Transfer available."""
set_available!(value::Transfer, val::Bool) = value.available = val
"""Set Transfer requirement."""
set_requirement!(value::Transfer, val::Float64) = value.requirement = val
"""Set Transfer ext."""
set_ext!(value::Transfer, val::Dict{String, Any}) = value.ext = val

InfrastructureSystems.set_forecasts!(value::Transfer, val::InfrastructureSystems.Forecasts) = value.forecasts = val
"""Set Transfer internal."""
set_internal!(value::Transfer, val::InfrastructureSystemsInternal) = value.internal = val


text:
abstract type TagBase end

struct DataSetConfig
    name::String
    basename::String
    datagen_backend_tag::TagBase
    io_backend_tag::TagBase
    nsteps::Int32
    step_factor::Int32
    compute_seconds::Float64
    streams::Vector{DataStreamConfig}
end

mutable struct DataSet
    cfg::DataSetConfig
    curr_step::Int32
    timestamp::Float64
    streams::Vector{DataStream}
end

DataSet(cfg::DataSetConfig) = DataSet(cfg, 1, time(), Float64[])


text:
export FDTD

struct FDTD

    ex::Array{Float64,2}
    ey::Array{Float64,2}
    bz::Array{Float64,2}

    function FDTD(mesh)

        nx, ny = mesh.nx, mesh.ny
        ex = zeros(nx, ny + 1)
        ey = zeros(nx + 1, ny)
        bz = zeros(nx, ny)

        new(ex, ey, bz)

    end

end

export faraday!

function faraday!(fdtd::FDTD, m::TwoDGrid, dt)

    nx, ny = m.nx, m.ny
    dx, dy = m.dx, m.dy

    for i = 1:nx, j = 1:ny
        dex_dy = (fdtd.ex[i, mod1(j + 1, ny)] - fdtd.ex[i, j]) / dy
        dey_dx = (fdtd.ey[mod1(i + 1, nx), j] - fdtd.ey[i, j]) / dx
        fdtd.bz[i, j] += dt * (dex_dy - dey_dx)
    end

end

export ampere_maxwell!

function ampere_maxwell!(fdtd::FDTD, m::TwoDGrid, jx, jy, dt)

    nx, ny = m.nx, m.ny
    dx, dy = m.dx, m.dy

    for i = 1:nx, j = 1:ny+1
        dbz_dy = (fdtd.bz[i, mod1(j, ny)] - fdtd.bz[i, mod1(j - 1, ny)]) / dy
        fdtd.ex[i, j] += dt * dbz_dy - dt * 0.5 * (jx[i, j] + jx[mod1(i + 1, nx), j])
    end

    for i = 1:nx+1, j = 1:ny
        dbz_dx = (fdtd.bz[mod1(i, nx), mod1(j, ny)] - fdtd.bz[mod1(i - 1, nx), j]) / dx
        fdtd.ey[i, j] -= dt * dbz_dx - dt * 0.5 * (jy[i, j] + jy[i, mod1(j + 1, ny)])
    end

end

export update_fields!

function update_fields!(ex, ey, bz, m::TwoDGrid, fdtd::FDTD)

    nx, ny = m.nx, m.ny

    for i = 1:nx+1, j = 1:ny+1
        ex[i, j] = 0.5 * (fdtd.ex[mod1(i - 1, nx), j] + fdtd.ex[mod1(i, nx), j])
        ey[i, j] = 0.5 * (fdtd.ey[i, mod1(j - 1, ny)] + fdtd.ey[i, mod1(j, ny)])
        bz[i, j] =
            0.25 * (
                fdtd.bz[mod1(i - 1, nx), mod1(j - 1, ny)] +
                fdtd.bz[mod1(i, nx), mod1(j - 1, ny)] +
                fdtd.bz[mod1(i - 1, nx), mod1(j, ny)] +
                fdtd.bz[mod1(i, nx), mod1(j, ny)]
            )
    end

end

export compute_energy

compute_energy(fdtd::FDTD, m::TwoDGrid) = 0.5 * log(sum(fdtd.ex .^ 2) * m.dx * m.dy)


text:
# Export Sections API
export Sections,
       getindex, length, lastindex, iterate, keys, eltype,
       findall, findfirst,
       handle, header, format_string

# Export Section API
export Section,
       deref, section_name, section_size, section_offset,
       section_address

# Export Datatypes
export SectionRef,
       read, seekstart, seek, eof, section_number

# Import Base methods for extension
import Base: read, seek, seekstart, eof, length, eltype, findall,
             findfirst, iterate, keys, lastindex

"""
    Sections

An abstraction over the concept of a collection of `Section` types within an
object file.  One can think of the `Sections` object containing the table of
section headers within the object file, whereas the `Section`/`SectionRef`
objects contain the actual section data itself.  The list of available API
operations is given below, with methods that subclasses must implement marked
in emphasis:

### Creation
  - *Sections()*

### Iteration
  - getindex()
  - *lastindex()*
  - length()
  - iterate()
  - keys()
  - eltype()

### Search
  - findall()
  - findfirst()

### Misc.
  - *handle()*
"""
abstract type Sections{H<:ObjectHandle} end

# Fairly simple iteration interface specification
@mustimplement lastindex(sections::Sections)
length(sections::Sections) = lastindex(sections)
keys(sections::Sections) = 1:length(sections)
iterate(sections::Sections, idx=1) = idx > length(sections) ? nothing : (sections[idx], idx+1)
eltype(::Type{S}) where {S <: Sections} = SectionRef

function getindex(sections::Sections, idx)
    # Punt off to `getindex_ref`
    oh = handle(sections)
    return getindex_ref(
        sections,
        section_header_offset(oh),
        section_header_size(oh),
        section_header_type(oh),
        SectionRef,
        idx
    )
end

"""
    findall(sections::Sections, name::String)

Return a list of sections that match the given `name`.
"""
function findall(sections::Sections, name::AbstractString)
    return findall(sections, [name])
end

"""
    findall(sections::Sections, name::String)

Return a list of sections that match one of the given `names`.
"""
function findall(sections::Sections, names::Vector{S}) where {S <: AbstractString}
    return [s for s in sections if section_name(s) in names]
end

"""
    findfirst(sections::Sections, name::String)

Return the first section that matches the given `name`.
"""
function findfirst(sections::Sections, name::AbstractString)
    return findfirst(sections, [name])
end

"""
    findfirst(sections::Sections, names::Vector{String})

Return the first section that matches on of the given `names`.
"""
function findfirst(sections::Sections, names::Vector{String})
    results = findall(sections, names)
    if isempty(results)
        error("Could not find any sections that match $(names)")
    end
    return first(results)
end

"""
    handle(sections::Sections)

Return the `ObjectHandle` that this `Sections` object belongs to
"""
@mustimplement handle(sections::Sections)

"""
    Section

An abstraction over the concept of a `Section` within an object file.  Because
many operations upon sections require global operations (access to the string
table, knowledge of position within the file, etc...) some operations are
defined only upon the `SectionRef` datatype.  As a user, the `SectionRef` type
should be the primary method of interacting with sections, as a developer
adding new object file formats, some methods must support `Section`s, others
must support only `SectionRef`s.  Note that any method that works on a
`Section` must also work with a `SectionRef`, see the `@derefmethod` macro for
a convenient helper macro to generate `SectionRef` -> `Section` wrapper
methods. The list of available API operations is given below, with methods that
subclasses must implement marked in emphasis:

### Creation:
  - *read()*

### Utility:
  - deref()

### IO-like operations:
  - read()

### Format-specific properties:
  - *section_name()*
  - *section_size()*
  - *section_offset()*
  - *section_address()*
"""
abstract type Section{H<:ObjectHandle} end

# read() is how we construct a `Section`
@mustimplement read(oh::ObjectHandle, ::Type{Section})

deref(section::Section) = section

"""
    read(oh::ObjectHandle, section::Section)

Read the contents of the section referred to by `section` from the given
`ObjectHandle`, returning a `Vector{UInt8}`.
"""
function read(oh::H, section::Section{H}) where {H<:ObjectHandle}
    # Seek to the section's location, then read it in!
    seek(oh, section_offset(section))
    return read(oh, section_size(section))
end

"""
    section_name(section::Section)

Return the name of the given section as a string.  In order to return a true
name, it is necessary to perform a lookup within the object's string table,
which cannot be done using just a `Section` object; use a `SectionRef` object
instead if you need that.  For sanity sake, this method will return a string,
but the contents of the string may be something like the offset within the
string table pointing to this `Section`'s name, e.g. "@strtab.123"
"""
@mustimplement section_name(section::Section)

"""
    section_size(section::Section)

The size of the actual data contained in the section. This should exclude any
padding mandated by the file format e.g. due to alignment rules
"""
@mustimplement section_size(section::Section)

"""
    section_offset(section::Section)

The offset of the section in the file, in bytes
"""
@mustimplement section_offset(section::Section)

"""
    section_address(section::Section)

The address of the section in virtual memory.
"""
@mustimplement section_address(section::Section)

"""
    SectionRef

Provides a reference to a `Section`, along with a reference to the
`ObjectHandle` this `Section` comes from.  This should be the primary method by
which users interact with sections inside object files.  The list of available
API operations is given below, with methods that subclasses must implement
marked in emphasis.  Note that this overlaps heavily with the `Section` object
API, this is by design as many of the methods are simply passthroughs to the
underlying `Section` API calls for ease of use.

### Creation:
  - *SectionRef()*

### Utility
  - *deref()*
  - *handle()*
  - *Sections()*

### IO-like operations:
  - read()
  - seekstart()
  - seek()
  - eof()

### Format-specific properties:
  - section_name()
  - *section_number()*
  - section_type()
  - section_size()
  - section_offset()
  - section_address()
"""
abstract type SectionRef{H<:ObjectHandle} end

"""
    SectionRef(sections::Sections, section::Section, idx)

Construct a `SectionRef` object pointing to the given `Section`, which itself
represents the `idx`'th section within the given `Sections`.
"""
@mustimplement SectionRef(sections::Sections, section::Section, idx)

"""
    deref(section::SectionRef)

Dereference the given `SectionRef` object to a `Section`.
"""
@mustimplement deref(section::SectionRef)

"""
    handle(section::SectionRef)

Return the `ObjectHandle` this `SectionRef` belongs to.  This method is
`SectionRef`-only.
"""
@mustimplement handle(section::SectionRef)

"""
    Sections(section::SectionRef)

Return the `Sections` collection this `section` belongs to.
"""
@mustimplement Sections(section::SectionRef)

"""
    section_number(section::SectionRef)

The index of the given section within the section header table.
"""
@mustimplement section_number(section::SectionRef)

@derefmethod section_name(section::SectionRef)
@derefmethod section_size(section::SectionRef)
@derefmethod section_offset(section::SectionRef)
@derefmethod section_address(section::SectionRef)

"""
    read(section::SectionRef)

Read the contents of the section referred to by `section`, returning a
`Vector{UInt8}`.
"""
function read(section::SectionRef)
    return read(handle(section), deref(section))
end

"""
    seekstart(section::SectionRef)

Seek to the beginning of `section` in the `ObjectHandle` it was loaded from.
"""
function seekstart(section::SectionRef)
    return seek(handle(section), section_offset(section))
end

"""
    seek(section::SectionRef, offset)

Seek to `offset` relative to `section` in the `ObjectHandle` that this
`SectionRef` refers to
"""
function seek(section::SectionRef, offset)
    return seek(handle(section), section_offset(section) + offset)
end

"""
    eof(section::SectionRef)

Returns `true` if the `ObjectHandle` that this `SectionRef` refers to has read
beyond the current section's extent
"""
function eof(section::SectionRef)
    # If the handle itself thinks it's an eof(), then yeah let's quit out
    if eof(handle(section))
        return true
    end

    # If we are beyond this section, then we also consider this the end of the
    # current section
    section_end = section_offset(section) + section_size(section)
    return position(handle(section)) >= section_end
end

# Common Printing
function show(io::IO, s::Union{Section{H},SectionRef{H}}) where {H <: ObjectHandle}
    print(io, "$(format_string(H)) Section")

    if !get(io, :compact, false)
        println(io)
        println(io, "       Name: $(section_name(s))")
        println(io, "       Size: 0x$(string(section_size(s), base=16))")
        println(io, "     Offset: 0x$(string(section_offset(s), base=16))")
        print(io,   "    Address: 0x$(string(section_address(s), base=16))")
    else
        print(io, " \"$(section_name(s))\"")
    end
end

show(io::IO, sects::Sections{H}) where {H <: ObjectHandle} = show_collection(io, sects, H)


text:
module PositiveFactorizations

using LinearAlgebra

export Positive

struct Positive{T<:Real} end

include("cholesky.jl")
include("eig.jl")

end # module


text:
## Issue examples ####

# Dimensions should be 330 x 570
cd("assets/")

wcpath = joinpath(ENV["SDMLAYERS_PATH"], "WorldClim", "BioClim", "10", "wc2.1_10m_bio_1.tif")
tmpfile = tempname()

query = `gdalwarp -te -145.0 20.0 -50.0 75.0 $(wcpath) $(tmpfile)`
run(query)

using ArchGDAL
d = ArchGDAL.read(tmpfile)
ArchGDAL.height(d) # 330
ArchGDAL.width(d) # 570

# hcat/vcat are wrong
l1 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=0.0, top=10.0);
l2 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=10.0, top=20.0);
l3 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=10.0, right=20.0, bottom=0.0, top=10.0);

vl1 = vcat(l1, l2) # latitudes span makes no sense
ml1 = hcat(l1, l3) # things look fine here, but are they?

latitudes(vl1)
vl1.bottom < vl1.top # false

latitudes(ml1)
ml1.bottom < ml1.top # false

# Inversed bounds
l1 = SimpleSDMPredictor(WorldClim, BioClim, 1)
l2 = SimpleSDMPredictor(copy(l1.grid), 180.0, -180.0, 90.0, -90.0) # looks the same, but bounds really are inversed

l1.grid == l2.grid # same grid
extrema(longitudes(l1)) == extrema(longitudes(l2)) # same longitude extremas (displayed by show method)
extrema(latitudes(l1)) == extrema(latitudes(l2)) # same latitude extremas too

longitudes(l1) == longitudes(l2) # false, inversed
latitudes(l1) == latitudes(l2) # false, inversed

# Geotiff writing
l1 = SimpleSDMPredictor(WorldClim, BioClim, 1)
length(l1)
geotiff(tempname(), l1)
length(l1)


text:
using PkgSRA
using Plots, Random, Flux, Turing, StatsPlots
pyplot()
Random.seed!(11)
include("paper_settings.jl");

#####
##### Define the controlled dynamical system
#####
# Load example problem
include(EXAMPLE_FOLDERNAME*"example_falling_ball.jl")
# Spatial forcing function
U_starts = [6.0]
    U_widths = 0.2;
    amplitude = 50.0;
U_func_kick(t, u) = U_func_time(t, u,
                                U_widths, U_starts,
                                F_dim=2,
                                amplitude=amplitude)
# Define the time dependent forcing function
U_func_wall(X) = U_func_spring(X, k=1e3, r=1.0);

# Solve the system
dyn_with_ctr = solve_ball_system(U_func_space=U_func_wall,
                                 U_func_time=U_func_kick)
# plot(dyn_with_ctr)
dat = Array(dyn_with_ctr)

# Get the true control signals
dyn_control_wall = zeros(size(dat))
for i in 1:size(dat, 2)
    dyn_control_wall[:,i] = U_func_wall(dat[:, i])
end

dyn_control_kick = zeros(size(dat))
for (i, t) in enumerate(ts)
    dyn_control_kick[:,i] = U_func_kick(t, dat[:, i])
end

# Take numerical gradient
numerical_grad = numerical_derivative(dat, ts)

#####
##### Bayesian residual process and control guess
#####

# Define model and sample
@model ball_residual(y, ind) = begin
    # Lorenz parameters
    g ~ Normal(10, 1.0)
    noise ~ Truncated(Normal(1, 1.0), 0, 20)

    t_vec = collect(ts)
    for i in train_ind
        du = ball_system(dat[:,i], g, t_vec[i])

        # y ~ du .+ MvNormal(length(du), noise)
        y[1,i] ~ Normal(du[1], noise)
        y[2,i] ~ Normal(du[2], noise)
    end
end;

# Better sampler: NUTS
iterations = 1000
    n_adapts = Int(iterations/5)
    j_max = 1.0
# Try to predict the GRADIENT from data
num_training_pts = 1000
    start_ind = 1
    train_ind = start_ind:num_training_pts+start_ind-1
    y = numerical_grad[:,train_ind]
chain = sample(ball_residual(y,train_ind), #noise, train_ind),
                NUTS(iterations, n_adapts, 0.6j_max));

plot(chain)

# Get a distribution of predictions
num_samples = 100
    param_samples = sample(chain, num_samples)
    vars = :g
    all_vals = zeros(num_samples, length(train_ind), size(dat,1))
    all_noise = zeros(num_samples)
for i in 1:num_samples
    these_params = param_samples[vars].value[i]
    all_noise[i] = param_samples[:noise].value[i]
    for (i_save, i_dat) in enumerate(train_ind)
        all_vals[i, i_save, :] = ball_system(dat[:,i_dat],
                                these_params, collect(ts)[i_dat])
    end
end

# Get the residuals
dat_grad = numerical_grad[:, train_ind]

# Calculate the residuals per variable
this_std = reshape(std(all_vals, dims=1), length(train_ind), size(all_vals,3))
    this_mean = reshape(mean(all_vals, dims=1), length(train_ind), size(all_vals,3))

# Final processing
# NOTE: This is the residual for the indices: train_ind
residual = dat_grad .- transpose(this_mean)
# ctr_guess = process_residual(residual, this_std[reconstruction_plot_ind])
ctr_guess = process_residual(residual, mean(all_noise))

plot(residual[2,:], ribbon=mean(all_noise))

plot(ctr_guess[2,:], label="Control Guess", lw=2)
    plot!(dyn_control_wall[2,:], label="Control Wall")
    plot!(dyn_control_kick[2,:], label="Control Kick")

#####
##### Fit first neural network: g(x)
#####
nn_dim = 16
    U_dim = 1 # Expect that the true dimensionality is 1
# Note: need to initialize the first layer so that the whole space is used
ctr_dyn = Chain(Dense(2,nn_dim, initb=(x)-> tspan[2].*rand(x)),
                Dense(nn_dim, nn_dim, σ),
                Dense(nn_dim, nn_dim, σ),
                Dense(nn_dim, U_dim),
                Dense(U_dim,2)) # Equivalent to the B matrix
ps = Flux.params(ctr_dyn)
    loss_U() = sum(abs2,ctr_dyn(dat) .- ctr_guess)
    loss_to_use = loss_U

# Initial fast learning
tol = 300
    dat_dummy = Iterators.repeated((), 100)
    opt = ADAM(1e-2)
    num_iter = 1
    max_iter = 5
while loss_to_use() > 5*tol
    Flux.train!(loss_to_use, ps, dat_dummy, opt, cb = ()->@show loss_to_use())
    global num_iter += 1
    num_iter > max_iter && break
end

plot(Flux.data(ctr_dyn(dat)[2,:]), label="g(x)", lw=2)
    plot!(dyn_control_wall[2,:], label="Control Wall")
    plot!(dyn_control_kick[2,:], label="Control Kick")

# plot(Flux.data(ctr_dyn(dat)[1,:]), label="g(x)", lw=2)
#     plot!(dyn_control_wall[1,:], label="Control Wall")
#     plot!(dyn_control_kick[1,:], label="Control Kick")

opt = ADAM(1e-5)
    num_iter = 1
while loss_to_use() > tol
    Flux.train!(loss_to_use, ps, dat_dummy, opt, cb = ()->@show loss_to_use())
    global num_iter += 1
    num_iter > max_iter && break
end
println("Finished learning control signal")

#####
##### Produce the plots
#####
# First: intrinsic dynamics
plot_data = plot(ts, dyn_with_ctr[1,:], label="Height", lw=3, legend=:bottomright);
    plot!(ts, dyn_with_ctr[2,:], label="Velocity", lw=3,
            legendfontsize=14, xticks=false, yticks=false);
    xlabel!("");
    title!("Observable Data", titlefontsize=24);

# Second: Controller
plot_control = plot(ts, vec(dyn_control_wall[2,:]), label="Ground",
                color=COLOR_DICT["control_true"], lw=3);
    plot_kick = plot!(ts, vec(dyn_control_kick[2,:]), label="Kick",
                color=COLOR_DICT["control_time"], lw=3,
                legendfontsize=14, xticks=false, yticks=false);
    xlabel!("Time", fontsize=16, legend=:bottomright);
    title!("Hidden Control Signals", titlefontsize=24);

# Creat the layout and plot
my_layout = @layout [p1; p2];
    p_final = plot(plot_data, plot_control, layout = my_layout)

# Save
fname = FIGURE_FOLDERNAME * "fig_spatial_and_temporal_control.png";
# savefig(p_final, fname)


text:
include("header.jl")

#=
# Manometry
In this notebook, we will learn
* How to use manometers to measure pressure
* How to combine several working fluids in a manometer
* How to calculate for tilted manometers
=#

# ### Set up the module
using MAE103
#-
using Plots

#=
In fluid dynamics, we are often interested in measuring pressure differences
between two points, often because those pressure differences allow us
to determine flow speed.

We'll discuss how we get flow speed from pressure difference later. For now, our
goal is to show how we can measure pressure difference between two
points. The key notion here is that we only need to measure pressure *difference*,
not absolute pressure. Also, one of those "points" might simply be the atmosphere,
where the pressure has some ambient value.

Manometers enable this by turning pressure difference into a measurable height of a
*working fluid* in the manometer. The manometer itself consists of a narrow tube,
bent into a U-shaped form and arranged vertically so that gravity and fluid statics can
be exploited. Sometimes the manometer has additional feature, e.g.,

- one of the channels in the tube might be at an angle rather than vertical
- one or more of the channels might be connected to a "cistern", with a
much wider diameter than the tube.

The working fluid is generally a **liquid**, and we can assume that it does not
change its volume during operation. This fact will be important in the analysis.
=#
#-
#=
### Example 1
We will demonstrate here how we use the fluid statics equation to determine pressure
by reading a measured height. Consider the example setup shown in the figure below:

#

We seek to measure the gauge pressure (i.e., relative to ambient) in the air in the tank
by measuring the heights of liquids. Suppose that the oil has specific gravity 0.9. Also,
mercury (Hg) has a very high specific gravity, 13.6: it is a liquid metal.

The heights shown are $h_1 = 36$ inches, $h_2 = 6$ inches, and $h_3 = 9$ inches.
Let's find the gauge pressure in the air, in psi.

First set the values:
=#
h1 = Height(36u"inch")
h2 = Height(6u"inch")
h3 = Height(9u"inch")
γoil = SpecificWeight(SpecificGravity(0.9))
γHg = SpecificWeight(SpecificGravity(13.6))
patm = Pressure(1u"atm");
#=
We work backward from the open end, where pressure is atmospheric. Just to be careful,
we will add this atmospheric pressure, so we are working with absolute pressures,
but then subtract it at the end to find our answer. At point 2, the pressure
increases by $γ_{\mathrm{Hg}}h_3$ from atmospheric:
=#
p2 = Pressure(patm + γHg*h3)
#=
Point 1 on the mercury side is at the same elevation as point 2, and in the same
fluid, so its pressures are equal.
=#
p1 = p2
#=
Also $p_1$ in oil is the same as in mercury. And this pressure is larger than
the pressure in the air by $\gamma_{\mathrm{oil}}(h_1 + h_2)$, so $p_{\mathrm{air}}$
=#
p_air = Pressure(p1 - γoil*(h1+h2))
#=
Remember, this was the absolute pressure. The gauge pressure is obtained by
subtracting the atmospheric pressure:
=#
p_airg = Pressure(p_air-patm)

#=
Finally, subtract the convert to psi
=#
value(p_airg,u"psi")

#=
In summary, the answer was

$$p_{\mathrm{air},g} = γ_{\mathrm{Hg}}h_3 - \gamma_{\mathrm{oil}}(h_1 + h_2)$$
=#

#=
In this example, the manometer used mercury. What if it had used water instead?
What would $h_3$ be (assuming the same gauge pressure in the air)? This is
given by

$$h_3 = \left(p_{\mathrm{air},g} + \gamma_{\mathrm{oil}}(h_1 + h_2)\right)/γ_w$$

where $γ_w$ is the specific weight of water. Let's find out what $h_3$ would be
(in feet).
=#
h3 = Height((p_airg + γoil*(h1+h2))/SpecificWeight(Water))
value(h3,u"ft")

#=
10.2 feet instead of 9 inches! There is big advantage in using a heavy working
liquid in a manometer, because it does not rise as high. It's easier to
get readings.
=#

#=
### Example 2
In the previous example, we pointed out that mercury is useful for translating
a large pressure into a relatively small height change. Sometimes we need the
opposite: if we want to measure a relatively small pressure difference,
we want it to be easy to read on the manometer. For this purpose, we can use an
inclined tube manometer, as in the setup shown in the figure below:

#

We read the length $l_2$ along the manometer tube. How does this length come into the
manometer equation? Only the vertical part of $l_2$ enters this equation, i.e.,
$l_2 \sin \theta$

The equation between $A$ and $B$ is

$$p_A + \gamma_1 h_1 - \gamma_2 l_2 \sin \theta - \gamma_3 h_3 = p_B$$

(From $A$, we drop down by $h_1$ with fluid 1, rise by $l_2 \sin \theta$ with fluid 2,
 then rise by $h_3$ with fluid 3.)

Let's assume that fluids 1 and 3 are gases, so their specific weights are negligible.
Then

$$p_B - p_A = \gamma_2 l_2 \sin \theta$$

or, in other words,

$$l_2 = \dfrac{p_B - p_A}{\gamma_2\sin\theta}$$

Let's try this out for a few different angles of inclined tube. We will assume that the
pressure difference is 100 Pa --- a pretty small pressure difference --- and the
working fluid (fluid 2) is water.
=#
dp = Pressure(100)
γ = SpecificWeight(Water)
#=
We'll start with a vertical tube (90 degrees)
=#
Θ = 90π/180
l = Length(dp/γ/sin(Θ))
value(l,u"cm")
#=
Now a tube at 45 degrees
=#
Θ = 45π/180
l = Length(dp/γ/sin(Θ))
value(l,u"cm")
#=
and finally a tube at 10 degrees
=#
Θ = 10π/180
l = Length(dp/γ/sin(Θ))
value(l,u"cm")
#=
Notice how the length goes up as the tube becomes more and more horizontal.
This makes it much easier to read the manometer without being prone to error.
=#


text:
#
println("Aq) Test of the Green(function) expansion.")

name          = "Lithium 1s^2 2s ground configuration"
refConfigs    = [Configuration("[He] 2s")]
greenSettings = GreenSettings(5, [0, 1, 2], 0.01, true, LevelSelection())
#
wa          = Representation(name, Nuclear.Model(8.), Radial.Grid(true), refConfigs,
                             ## GreenExpansion( Atomic.SingleCSFwithoutCI(), Basics.DeExciteSingleElectron(),
                             ## GreenExpansion( Atomic.CoreSpaceCI(), Basics.DeExciteSingleElectron(),
                                GreenExpansion( AtomicState.DampedSpaceCI(), Basics.DeExciteSingleElectron(),
                                             [LevelSymmetry(1//2, Basics.plus), LevelSymmetry(3//2, Basics.plus)], 3, greenSettings) )
println(wa)

wb = generate(wa, output=true)


text:
#=
Copyright 2018 INSIGNEO Institute for in silico Medicine

Licensed 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.
=#

"""
    checkConvergence(edge_list, vessels :: Array{Vessel, 1})

Compute the maximum error in the pressure and volumetric flow rate waveforms between two
cardiac cycles at the midpoint of all the vessels in the network.
"""
function checkConvergence(edge_list, vessels :: Array{Vessel, 1})
    err = zeros(size(edge_list)[1],2) .+ 100

    @inbounds for i in 1:size(edge_list)[1]
        v = vessels[i]
        lbl = v.label

        w_last = v.Q_l
        w_temp = v.Q_t

        err[i,1] = maximum(abs.((w_last[:,4].-w_temp[:,4])./w_last[:,4])*100)

        w_last = v.P_l
        w_temp = v.P_t

        err[i,2] = maximum(abs.((w_last[:,4].-w_temp[:,4])./w_last[:,4])*100)
    end

    return maximum(err)
end


text:
abstract type AbstractPattern end

include("patterns/empirical_pattern.jl")


text:
module MatUtils

import Base: copy, real, imag, map, size
import LinearAlgebra: diag

export zeroOut
export diag, diagProduct, minor, hinv

function zeroOut(x::Real,ϵ=1e-14)
    x < ϵ ? x = 0 : x
end

function zeroOut(x::Complex,ϵ=1e-14)
    xr = real(x); xi = imag(x)
    return zeroOut(xr,ϵ) + im*zeroOut(xi,ϵ)
end

function zeroOut(x::AbstractArray,ϵ=1e-14)
    map(x->zeroOut(x,ϵ),x)
end

function diag(A,offset::Integer)
    n,m = size(A)
    d = diag(A)

    if offset < 0
        offset = -offset
        @assert offset < n
        d = diag(A[1+offset:8,1:m])
    end

    if offset > 0
        @assert offset < m
        d = diag(A[1:n,1+offset])
    end

    return d
end

#gives the product of a diagonal of matrix A
#from index a (top left) to b (bottom right)
function diagProduct(A,offset,a,b)
    return reduce(*,diag(A,offset)[a:b])
end

#returns a matrix with row "row" and
#column "col" deleted.  A zero for row
#or column implies no deleted row or
#column respectively.
function minor(A,row=0,column=0)
    n,m = size(A)
    B = A[1:n .!= row,:] #remove row
    C = B[:,1:m .!= column] #remove col
    return C
end

hinv(i,j) = ((-1)^(i+j)*det(minor(H,j,i)))/det(H)

end #module MatUtils

max_stars_repo_name: sphanna/MatUtils

max_issues_repo_name: sphanna/MatUtils

max_forks_repo_name: sphanna/MatUtils


text:
using Documenter
using Registrator

makedocs(modules=[Registrator],
         sitename="Registrator.jl",
         pages=["Home" => "index.md",
                "Hosting Your Own" => "hosting.md",
                "Using Docker" => "docker.md",
                "Comment Bot" => "commentbot.md",
                "Web UI" => "webui.md"])

deploydocs(repo="github.com/JuliaRegistries/Registrator.jl.git")


text:
"""
Placeholder for a short summary about NonhomotheticCES.
"""
module NonhomotheticCES

export NonhomotheticCESUtility, log_consumption_aggregator, log_sectoral_consumptions

using ArgCheck: @argcheck
using DocStringExtensions: SIGNATURES
using LinearAlgebra: dot
using LogExpFunctions: logsumexp
using Requires: @require
using StaticArrays: SVector
using Statistics: mean
using UnPack: @unpack

include("internals.jl")

struct NonhomotheticCESUtility{N,Tσ<:Real,TΩ̂<:Real,Tϵ<:Real}
    σ::Tσ
    Ω̂s::SVector{N,TΩ̂}
    ϵs::SVector{N,Tϵ}
    function NonhomotheticCESUtility(σ::Tσ, Ω̂s::SVector{N,TΩ̂},
                   ϵs::SVector{N,Tϵ}) where {N,Tσ<:Real,TΩ̂<:Real,Tϵ<:Real}
        check_σ_ϵs(σ, ϵs)
        new{N,Tσ,TΩ̂,Tϵ}(σ, Ω̂s, ϵs)
    end
end

Broadcast.broadcastable(U::NonhomotheticCESUtility) = Ref(U)

"""
$(SIGNATURES)

Non-homothetic CES preferences, as defined in

*Comin, D., Lashkari, D., & Mestieri, Martí (2021). Structural change with long-run
 income and price effects. Econometrica, 89(1), 311–374.*

### Definition

Let ``E = ∑ᵢ pᵢ Cᵢ`` denote that *total expenditure*. The consumption aggregator `C` is
defined implicitly by

```math
E^{1 - \\sigma} = \\sum_i \\Omega_i (C^{\\epsilon_i} p_i)^{1 - \\sigma}
```

In the actual calculation and parametrization, we use logs (``Ê = log(E)`` etc) for improved
floating point accuracy.

### Arguments

- `σ`: elasticity of substitution between goods of different sectors, `> 0`, `≠ 1`.
- `Ω̂s`: **log** of sector weights
- `ϵs`: sectoral non-homotheticity parameters.

### Suggestions

- Use `SVector`s for vector arguments whenever possible.

- Normalize to `ϵ = 1` and `Ω̂ = 0` for a *base good* as described in equation (10) of the
  paper. This is not mandatory, but improves numerical performance.
"""
function NonhomotheticCESUtility(σ::Real, Ω̂s::AbstractVector, ϵs::AbstractVector)
    N = length(Ω̂s)
    @argcheck N == length(ϵs)
    NonhomotheticCESUtility(σ, SVector{N}(Ω̂s), SVector{N}(ϵs))
end
# FIXME: SizedVector constructors

"""
$(SIGNATURES)

Calculate the **log** consumption aggregator for the given utility, with **log** sector
prices `p̂s` and **log** expenditure `Ê`; within (absolute) tolerance `tol`.
"""
function log_consumption_aggregator(U::NonhomotheticCESUtility{N,Tσ,TΩ̂,Tϵ},
                                    p̂s::SVector{N,Tp̂}, Ê::TÊ;
                                    tol = DEFAULT_TOL) where {N,Tσ,TΩ̂,Tϵ,TÊ,Tp̂}
    @unpack σ, Ω̂s, ϵs = U
    @argcheck isfinite(Ê) && all(isfinite, p̂s) DomainError
    calculate_Ĉ(Ê, σ, Ω̂s, ϵs, p̂s; tol = tol, skipcheck = true)
end

function __init__()
    @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" begin
        using .ForwardDiff: Dual, value, partials, valtype
        import NonhomotheticCES: calculate_Ĉ

        partials_product(x, y, α = 1) = mapreduce((x, y) -> partials(x) * y * α, +, x, y)

        function calculate_Ĉ(::Type{T}, Ê, σ, Ω̂s, ϵs, p̂s;
                             tol = DEFAULT_TOL, skipcheck = false) where {T<:Dual}
            vÊ, vσ, vΩ̂s, vϵs, vp̂s = value(Ê), value(σ), value.(Ω̂s), value.(ϵs), value.(p̂s)
            Ĉ = calculate_Ĉ(valtype(T), vÊ, vσ, vΩ̂s, vϵs, vp̂s;
                            tol = tol, skipcheck = skipcheck)
            Zs, ∑Zϵ = calculate_Zs_∑Zϵ(Ĉ, vÊ, vσ, vΩ̂s, vϵs, vp̂s)
            # NOTE: we could be clever here and only calculate what is actually needed
            ∂p̂s = calculate_∂p̂s(Zs, ∑Zϵ)
            dĈ = (partials(Ê) * calculate_∂Ê(Zs, ∑Zϵ) +                      # ∂Ê
                  partials(σ) * calculate_∂σ(Zs, ∑Zϵ, Ĉ, vÊ, vσ, vϵs, vp̂s) + # ∂σ
                  partials_product(Ω̂s, ∂p̂s, inv(1 - vσ)) +                   # ∂Ω̂
                  partials_product(ϵs, ∂p̂s, Ĉ) +                             # ∂ϵs
                  partials_product(p̂s, ∂p̂s))                                 # ∂p̂s
            T(Ĉ, dĈ)
        end
    end
end

"""
$(SIGNATURES)

Calculate **log** sectoral consumptions, return as a vector. Arguments are in logs; see
[`log_consumption_aggregator`](@ref) which whould also be used to obtain `Ĉ`.

The budget constraint holds, ie
```julia
logsumexp(log_sectoral_consumptions(U, p̂s, Ê, Ĉ) .+ p̂s) ≈ Ê
```
"""
function log_sectoral_consumptions(U::NonhomotheticCESUtility, p̂s, Ê, Ĉ)
    @unpack σ, Ω̂s, ϵs = U
    Ω̂s .- σ .* (p̂s .- Ê) + ((1 - σ) * Ĉ) .* ϵs
end

end # module


text:
using PowerSystems
using JuMP
using PowerSimulations
using GLPK
const PS = PowerSimulations

# ED Testing
base_dir = dirname(dirname(pathof(PowerSystems)))
include(joinpath(base_dir,"data/data_5bus.jl"))
sys5 = PowerSystem(nodes5, generators5, loads5_DA, branches5, nothing, 100.0,runchecks=false);
simple_reserve = PowerSystems.StaticReserve("test_reserve",sys5.generators.thermal,60.0,[gen.tech for gen in sys5.generators.thermal])

# ED with thermal gen, static load, copper plate
@test try
    @info "ED with thermal gen, static load, copper plate"
    ED = PS.PowerOperationModel(PS.EconomicDispatch,
                            [(device = ThermalGen, formulation =PS.ThermalDispatch)],
                            nothing,
                            [(device=Line, formulation=PS.PiLine)],
                            PS.CopperPlatePowerModel,
                            nothing,
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(ED,sys5)
    JuMP.optimize!(ED.model,with_optimizer(GLPK.Optimizer))
    (ED.model.moi_backend.model.optimizer.termination_status == JuMP.MOI.Success)  ? true : @error("solver returned with nonzero status")
true finally end

# ED with thermal and curtailable renewable gen, static load, copper plate
@test try
    @info "ED with thermal and curtailable renewable gen, static load, copper plate"
    ED = PS.PowerOperationModel(PS.EconomicDispatch,
                            [(device = ThermalGen, formulation =PS.ThermalDispatch),
                             (device = RenewableGen, formulation = PS.RenewableCurtail)],
                            nothing,
                            nothing,
                            [(device=Line, formulation=PS.PiLine)],
                            PS.CopperPlatePowerModel,
                            nothing,
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(ED,sys5)
    JuMP.optimize!(ED.model,with_optimizer(GLPK.Optimizer))
    (ED.model.moi_backend.model.optimizer.termination_status == JuMP.MOI.Success) ? true : @error("solver returned with nonzero status")
true finally end

# ED with thermal and fixed renewable gen, interruptable load, copper plate
@test try
    @info "ED with thermal and fixed renewable gen, interruptable load, copper plate"
    ED = PS.PowerOperationModel(PS.EconomicDispatch,
                            [(device = ThermalGen, formulation = PS.ThermalDispatch),
                             (device = RenewableGen, formulation = PS.RenewableCurtail)],
                            [(device = ElectricLoad, formulation = PS.InterruptibleLoad)],
                            nothing,
                            [(device=Line, formulation=PS.PiLine)],
                            PS.CopperPlatePowerModel,
                            nothing,
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(ED,sys5)
    JuMP.optimize!(ED.model,with_optimizer(GLPK.Optimizer))
    (ED.model.moi_backend.model.optimizer.termination_status == JuMP.MOI.Success)  ? true : @error("solver returned with nonzero status")
true finally end

# ED with thermal gen, copper plate, and reserve
@test try
    @info "ED with thermal gen, copper plate, and reserve"
    ED = PS.PowerOperationModel(PS.EconomicDispatch,
                            [(device = ThermalGen, formulation =PS.ThermalDispatch)],
                            nothing,
                            [(device=Line, formulation=PS.PiLine)],
                            PS.CopperPlatePowerModel,
                            [(service = simple_reserve, formulation = PS.RampLimitedReserve)],
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(ED,sys5)
    JuMP.optimize!(ED.model,with_optimizer(GLPK.Optimizer))
    (ED.model.moi_backend.model.optimizer.termination_status == JuMP.MOI.Success)  ? true : @error("solver returned with nonzero status")
true finally end

# ED with thermal gen, PTDF
@test try
    @info "ED with thermal gen, copper plate, and reserve"
    ED = PS.PowerOperationModel(PS.EconomicDispatch,
                            [(device = ThermalGen, formulation =PS.ThermalDispatch),
                            (device = RenewableGen, formulation = PS.RenewableCurtail)],
                            nothing,
                            nothing,
                            [(device=Branch, formulation=PS.PiLine)],
                            PS.StandardPTDF,
                            nothing,
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(sys5,ED)
    JuMP.optimize!(ED.model,with_optimizer(GLPK.Optimizer))
    (ED.model.moi_backend.model.optimizer.termination_status == JuMP.MOI.Success)  ? true : @error("solver returned with nonzero status")
true finally end

# UC Testing
base_dir = dirname(dirname(pathof(PowerSystems)))
include(joinpath(base_dir,"data/data_5bus_uc.jl"))
sys5 = PowerSystem(nodes5, generators5, loads5_DA, branches5, nothing, 100.0, runchecks=false);
simple_reserve = PowerSystems.StaticReserve("test_reserve",sys5.generators.thermal,60.0,[sys5.generators.thermal[1].tech])

# UC with thermal gen, static load, copper plate
@test try
    @info "UC with thermal gen, static load, copper plate"
    UC = PS.PowerOperationModel(PS.UnitCommitment,
                            [(device = ThermalGen, formulation =PS.StandardThermalCommitment)],
                            nothing,
                            nothing,
                            [(device=Branch, formulation=PS.PiLine)],
                            PS.CopperPlatePowerModel,
                            nothing,
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(UC,sys5)
    JuMP.optimize!(UC.model,with_optimizer(GLPK.Optimizer))
    (UC.model.moi_backend.model.optimizer.termination_status == JuMP.MOI.Success)  ? true : @error("solver returned with nonzero status")
true finally end

# UC with thermal and curtailable renewable gen, static load, copper plate
@test try
    @info "UC with thermal and curtailable renewable gen, static load, copper plate"
    UC = PS.PowerOperationModel(PS.EconomicDispatch,
                            [(device = ThermalGen, formulation =PS.StandardThermalCommitment),
                             (device = RenewableGen, formulation = PS.RenewableCurtail)],
                            nothing,
                            nothing,
                            [(device=Line, formulation=PS.PiLine)],
                            PS.CopperPlatePowerModel,
                            nothing,
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(UC,sys5)
    JuMP.optimize!(UC.model,with_optimizer(GLPK.Optimizer))
true finally end

# UC with thermal and fixUC renewable gen, interruptable load, copper plate
@test try
    @info "UC with thermal and fixUC renewable gen, interruptable load, copper plate"
    UC = PS.PowerOperationModel(PS.UnitCommitment,
                            [(device = ThermalGen, formulation = PS.StandardThermalCommitment),
                             (device = RenewableGen, formulation = PS.RenewableCurtail)],
                            [(device = ElectricLoad, formulation = PS.InterruptibleLoad)],
                            nothing,
                            [(device=Line, formulation=PS.PiLine)],
                            PS.CopperPlatePowerModel,
                            nothing,
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(UC,sys5)
    JuMP.optimize!(UC.model,with_optimizer(GLPK.Optimizer))
    (UC.model.moi_backend.model.optimizer.termination_status == JuMP.MOI.Success)  ? true : @error("solver returned with nonzero status")
true finally end

# UC with thermal gen, copper plate, and reserve
@test try
    @info "UC with thermal gen, copper plate, and reserve"
    UC = PS.PowerOperationModel(PS.EconomicDispatch,
                            [(device = ThermalGen, formulation =PS.StandardThermalCommitment),
                            (device = RenewableGen, formulation = PS.RenewableCurtail)],
                            nothing,
                            [(device=Line, formulation=PS.PiLine)],
                            PS.CopperPlatePowerModel,
                            [(service = simple_reserve, formulation = PS.RampLimitedReserve)],
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(sys5,UC)
    JuMP.optimize!(UC.model,with_optimizer(GLPK.Optimizer))
    (UC.model.moi_backend.model.optimizer.termination_status == JuMP.MOI.Success)  ? true : @error("solver returned with nonzero status")
true finally end

# UC with thermal gen, copper plate, and PTDF
@test try
    @info "UC with thermal gen, copper plate, and reserve"
    UC = PS.PowerOperationModel(PS.EconomicDispatch,
                            [(device = ThermalGen, formulation =PS.StandardThermalCommitment),
                            (device = RenewableGen, formulation = PS.RenewableCurtail)],
                            nothing,
                            nothing,
                            [(device=Branch, formulation=PS.PiLine)],
                            PS.StandardPTDF,
                            nothing,
                            sys5,
                            Model(),
                            false,
                            nothing)
    PS.buildmodel!(UC,sys5)
    JuMP.optimize!(UC.model,with_optimizer(GLPK.Optimizer))
    (UC.model.moi_backend.model.optimizer.termination_status == JuMP.MOI.Success)  ? true : @error("solver returned with nonzero status")
true finally end


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

sojurn_times = Float64[]

number_in_system_integral = 0
prev_number_in_system = 0
prev_time = 0

processed_sojurn_times = []

function mode_1(time::Float64, state::NetworkState, event::TimedEvent)
    global number_in_system_integral += (time - prev_time) * prev_number_in_system
    global prev_number_in_system = state.number_in_system
    global prev_time = time
    return nothing
end

function mode_2(time::Float64, state::NetworkState, event::TimedEvent)
    state.number_in_system_decreased && push!(sojurn_times, time - event.arrival_time)
    return nothing
end

modes_dict = Dict(1 => mode_1, 2 => mode_2)

function run_script(params::NetworkParameters; max_time::Float64, lambda_range::Vector{Float64}, mode::Int=1, plot_num::Int=1)

    global modes_dict

    plot_data_callbacks = modes_dict[mode]

    mean_system_job_totals = zeros(Float64, length(lambda_range))
    proportion_orbiting_jobs = zeros(Float64, length(lambda_range))

    for i in 1:length(lambda_range)

        global sojurn_times, prev_time, prev_number_in_system, number_in_system_integral, processed_sojurn_times

        params.λ = lambda_range[i]

        Random.seed!(1)

        println("Starting simulation for λ = ", lambda_range[i])

        sim_states = simulate(NetworkState(0, zeros(Int8, params.L), 0), TimedEvent(ArrivalEvent(),0.0, 0, 0.0), params, max_time = max_time, callback = plot_data_callbacks)

        println("Simulation for lambda = ", lambda_range[i], " complete. Processing data for lambda = ", lambda_range[i], " now.")

        if mode == 1
            # Logic to generate the data needed for Plot 1.
            mean_system_job_totals[i] = number_in_system_integral/max_time

            # Logic to generate the data needed for Plot 2.
            proportion_data = sim_states.orbiting_jobs/sim_states.number_in_system
            !isnan(proportion_data) ? proportion_orbiting_jobs[i] = proportion_data : proportion_orbiting_jobs[i] = 0.0

            # proportion_orbiting_jobs[i] = sim_states.orbiting_jobs/sim_states.number_in_system == NaN ? 0.0 : sim_states.orbiting_jobs/sim_states.number_in_system

        elseif mode == 2
            # Logic to generate the data needed for Plot 3.
            push!(processed_sojurn_times, sojurn_times)
        else
            throw(DomainError(mode, "Mode number currently unsupported."))
        end

        println("Data for lambda = ", lambda_range[i], " has now been processed.\n")

        sojurn_times = Float64[]
        prev_time = 0
        prev_number_in_system = 0
        number_in_system_integral = 0
    end
    scenario_plots(processed_sojurn_times, mean_system_job_totals, proportion_orbiting_jobs, plot_num, lambda_range, mode)
end

#################################
### NetworkParameters Outline ###
#################################

"""
@with_kw mutable struct NetworkParameters
    L::Int #amount of stations
    gamma_shape::Float64 #external rate of arrivals
    λ::Float64 #This is undefined for the scenarios since it is varied
    η::Float64 #overflow movement rate
    μ_vector::Vector{Float64} #service rates
    P::Matrix{Float64} #routing matrix
    Q::Matrix{Float64} #overflow matrix
    p_e::AbstractWeights #external arrival distribution
    K::Vector{Int} #-1 means infinity
end
"""

###############################
### HOW TO USE THIS PROGRAM ###
###############################

"""
Please follow the instructions to README.md which explains how to use this program. The following variables test_params and lambdas along with the subsequent call to run_script are all the inputs needed to run this program. Once the inputs are what you want to them to be, run this script without debugging in VSCode with Ctrl+F5 and wait for the script to finish processing the request. Once finished, the script will output PNG files related to the mode that you have selected.
"""

test_params = NetworkParameters(  L=3,
                                gamma_shape = 3.0,
                                λ = NaN,
                                η = 4.0,
                                μ_vector = ones(3),
                                P = [0 1.0 0;
                                    0 0 1.0;
                                    0.5 0 0],
                                Q = [0 0.5 0;
                                     0 0 0.5;
                                     0.5 0 0],
                                p_e = AnalyticWeights([1.0, 0, 0]),
                                K = fill(5,3))

lambdas = [i for i in 0.2:0.2:1.2]

run_script(test_params, max_time = 10000.0, lambda_range = lambdas, mode=2, plot_num=3)


text:
module MusicAlbums

export Model, Mapper, Services, Resources, Client

include("connectionPools.jl")
using .ConnectionPools

include("workers.jl")
using .Workers

include("model.jl")
using .Model

include("auth.jl")
using .Auth

include("contexts.jl")
using .Contexts

include("mapper.jl")
using .Mapper

include("services.jl")
using .Service

include("resources.jl")
using .Resource

include("../Client/client.jl")
using .Client

function run(dbfile, authkeysfile)
    Workers.init()
    Mapper.init(dbfile)
    Auth.init(authkeysfile)
    Resource.run()
end

end # module

max_stars_repo_name: vanhp/MicroServiceJulia

max_issues_repo_name: vanhp/MicroServiceJulia

max_forks_repo_name: vanhp/MicroServiceJulia


text:
function isequal(S::AbstractStimulus, Q::AbstractStimulus)
    typeof(S) == typeof(Q) || return false
    for (fp,fq) in zip(fieldnames(typeof(S)),fieldnames(typeof(Q)))
        isequal(getfield(S, fp), getfield(Q, fq)) || return false
    end
    return true
end

function hash(S::AbstractStimulus, h::UInt)
    _h = h
    for fn in fieldnames(S)
        _h = hash(getfield(S, fn), _h)
    end
    return _h
end

################################################################################
#### File IO
################################################################################
"""
    savestimulus(S::AbstractStimulus; [filename], [dir], [old])

Save the given stimulus to disk using the `JLD` package. Returns the full path to the saved
file. Deletes file `old` if it exists in `dir`. By default `old` is `nothing` to avoid
unnecessarily deleting and rewriting files.

Default `filename` is `string(hash(S))`. Default `dir` is `Pkg.dir(GrayScaleStimuli)/saved`.

If `dir` doesn't exist, will use `mkpath`. If file exists, contents will be overwritten.
"""
function savestimulus(S::AbstractStimulus;
    filename=string(hash(S)),
    dir=joinpath(Pkg.dir("GrayScaleStimuli"), "saved"),
    old=nothing)

    # _dir = abspath(dir)
    if !ispath(dir)
        mkpath(dir)
    end
    if old != nothing
        _old = endswith(old, ".jld") ? old : old * ".jld"
        if isfile(joinpath(dir, _old))
            rm(joinpath(dir, _old))
        end
    end
    _fn = endswith(filename, ".jld") ? filename : filename * ".jld"
    # save(joinpath(_dir, _fn), "S", S)
    jldopen(joinpath(dir, _fn), "w") do file
        addrequire(file, GrayScaleStimuli)
        write(file, "S", S)
    end
    return joinpath(dir, _fn)
end

"""
    loadstimulus(filename; [dir])
    loadstimulus(h::UInt; [dir])

Returns the stimulus saved in `filename` (or identified by hash `h`) as saved by
`savestimulus`. Default `dir` is `Pkg.dir(GrayScaleStimuli)/saved`

"""
function loadstimulus(filename; dir=joinpath(Pkg.dir("GrayScaleStimuli"), "saved"))
    _fn = endswith(filename, ".jld") ? filename : filename * ".jld"
    _fn = basename(_fn)
    return load(joinpath(dir, _fn), "S")
end
loadstimulus(h::UInt; dir=joinpath(Pkg.dir("GrayScaleStimuli"),"saved")) = loadstimulus(string(h); dir=dir)


text:
#From http://home.online.no/~pjacklam/notes/invnorm/

function FMΦInv(p::Float64)::Float64

    a1::Float64 = -39.6968302866538
    a2::Float64 = 220.946098424521
    a3::Float64 = -275.928510446969
    a4::Float64 = 138.357751867269
    a5::Float64 = -30.6647980661472
    a6::Float64 = 2.50662827745924

    b1::Float64 = -54.4760987982241
    b2::Float64 = 161.585836858041
    b3::Float64 = -155.698979859887
    b4::Float64 = 66.8013118877197
    b5::Float64 = -13.2806815528857

    c1::Float64 = -7.78489400243029E-03
    c2::Float64 = -0.322396458041136
    c3::Float64 = -2.40075827716184
    c4::Float64 = -2.54973253934373
    c5::Float64 = 4.37466414146497
    c6::Float64 = 2.93816398269878

    d1::Float64 = 7.78469570904146E-03
    d2::Float64 = 0.32246712907004
    d3::Float64 = 2.445134137143
    d4::Float64 = 3.75440866190742

    #Rational approximation for lower region.

    if p < 0.02425
        q::Float64 = √(-2.0 * log(p))
        ans::Float64 = (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1.0)
    elseif p < 0.97575
        q = (p - 0.5) * (p - 0.5)
        ans = (((((a1 * q + a2) * q + a3) * q + a4) * q + a5) * q + a6) * (p - 0.5) / (((((b1 * q + b2) * q + b3) * q + b4) * q + b5) * q + 1)
    else
        q = √(-2 * log(1 - p))
        ans = -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1)
    end

    return ans
end

FMΦ(z::Real)::Real = 0.5+0.5*erf(z*sqrt2Inv)
FMϕ(z::Real)::Real = sqrt2PiInv*exp(-z^2.0/2.0)

#vcats like with rbind in r
function vbind(df1::AbstractDataFrame, df2::AbstractDataFrame)::AbstractDataFrame
  #get info on what we are joining
  fields1::Vector{<:DField} = names(df1)
  fields2::Vector{<:DField} = names(df2)
  allFields::Vector{<:DField} = union(fields1, fields2)

  rows1::Int = size(df1,1)
  rows2::Int = size(df2,1)

  for f::DField ∈ setdiff(allFields, fields1)
    origType::Type = eltype(df2[f])

    if !(Missing <: origType) #in this case we need to promote the eltype
      origType = Union{origType, Missing}
      df2[f] = Vector{origType}(df2[f])
    end

    df1[ f] = Vector{origType}(missing, rows1)
  end

  for f::DField ∈ setdiff(allFields, fields2)
    origType = eltype(df1[f])

    if !(Missing <: origType) #in this case we need to promote the eltype
      origType = Union{origType, Missing}
      df1[f] = Vector{origType}(df1[f])
    end

    #println("origType: $origType")
    df2[ f] = Vector{origType}(missing, rows2)
  end

  return [df1; df2]
end

#recursive version that works similar to vcat with the above stipulations
function vbind(dfs::T...)::T where T<:AbstractDataFrame
  if length(dfs) == 1
    return dfs[1]
  else
    pivot::Int = length(dfs) ÷ 2
    return vbind(vbind(dfs[1:pivot]...), vbind(dfs[(pivot+1):end]...))
  end
end

max_stars_repo_name: clintonTE/Finometrics

max_issues_repo_name: clintonTE/Finometrics

max_forks_repo_name: clintonTE/Finometrics


text:
"""
    radar(names, values)

Creates an `EChart` where region inside plotted line optionally filled with color.

## Methods
```julia
radar(names::AbstractVector, values::AbstractVector{<:Union{Missing, Real}},
    max::AbstractVector{<:Union{Missing, Real}})
radar(names::AbstractVector, values::AbstractArray{<:Union{Missing, Real},2},
	max::AbstractVector{<:Union{Missing, Int, AbstractFloat, Rational}})
```

## Arguments
* `fill::Union{Bool, AbstractVector} = false` : fill area inside marks with color?
* `legend::Bool` : display legend?
* `kwargs` : varargs to set any field of resulting `EChart` struct

## Notes

## Examples
```julia
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
ar = area(x, y)
```
"""
function radar(names::AbstractVector,
			values::AbstractVector{<:Union{Missing, Real}},
			max::AbstractVector{<:Union{Missing, Real}};
            fill::Bool = false,
			legend::Bool = false,
			kwargs...)

	#Validate arrays are same length
	if size(names)[1] != size(values)[1]
		error("Arrays names and values need to have the same length.")
	end

	ec = newplot(kwargs, ec_charttype = "radar")

	ec.radar = [Radar(indicator = arrayofdicts(name = names, max = max))]

	ec.series = [RadarSeries(name = "Series 1", data = [Dict{Any, Any}("value" => values)])]

	# Fill area if requested
	fill!(ec, 1, fill)

	#Add legend if requested
	legend ? legend!(ec) : nothing

	return ec

end

function radar(names::AbstractVector,
			   values::AbstractArray{<:Union{Missing, Real},2},
			   max::AbstractVector{<:Union{Missing, Real}};
	 			fill::Union{Bool, AbstractVector} = false,
				legend::Bool = true,
				kwargs...)

	# Call 1-D method to build base
	ec = radar(names, values[:,1], max; kwargs...)

	# Append remaining Y data
	for i in 2:size(values)[2]
		push!(ec.series, RadarSeries(data = [Dict{Any, Any}("value" => values[:,i])]))
	end

	# Fill area if requested
	# This is smelly
	typeof(fill) == Bool ? fill = [fill for i in 1:size(values)[2]] : fill = fill
	cols = length(fill)
	fill!(ec, cols, fill)

	seriesnames!(ec)

	#Add legend if requested
	legend ? legend!(ec) : nothing

	return ec

end


text:
# this is a sorting algorithm where the thing to be sorted and an index is sorted
# by packing the two into a 64bit word
using Revise
using SortingAlgorithms, SortingLab
using InternedStrings

@time samplespace = InternedString.("id".*dec.(1:1_000_000, 10));
srand(1);
@time svec = rand(samplespace, 100_000_000);

# @time isvec = InternedString.(svec)

#@time SortingAlgorithms.sortperm(svec, alg=StringRadixSort); # 51-55
# @time sortperm(svec); # 228
# @time unsafe_load.(Ptr{UInt}.(pointer.(svec)) .+ 4);
# @time ptrx = Ptr{UInt}.(pointer.(svec));
# @time unsafe_load.(ptrx .+ 4);
@time fsortperm(["abc","def"]);

srand(1);
svec = rand(samplespace, 100_000_000);
fsortperm_timing = @elapsed ss = fsortperm(svec);
issorted(svec[ss])

srand(1);
svec = rand(samplespace, 100_000_000);
sa_timing = @elapsed ss = SortingAlgorithms.sortperm_radixsort(svec);
issorted(svec[ss])

using RCall

r_timing = R"""
ss = sprintf("id%010d",1:1e6)
x = sample(ss, 1e8, replace=T)
system.time(order(x, method="radix"))
"""[3]

srand(1);
svec = rand(samplespace, 100_000_000);
base_timing = @elapsed ss = SortingAlgorithms.sortperm(svec);
issorted(svec[ss])

using StatPlots

secs_bar = bar(["fsortperm", "SortingAlgorithms.jl (unmerged PR)", "Julia baase", "R radixsort"],
    [fsortperm_timing, sa_timing, base_timing, r_timing])

rela_bar = bar(["fsortperm", "SortingAlgorithms.jl (unmerged PR)", "Julia base", "R radixsort"],
    [fsortperm_timing/r_timing, sa_timing/r_timing, base_timing/r_timing, r_timing/r_timing],
    label="multiple of R timing", title="Sortperm (R's order) performance",
    yticks = [0,5,10,15,20,25,30,35,40])

savefig(rela_bar,"sortperm_perf.png")

plot(secs_bar, rela_bar)