text:
struct DESPOT{S,A,O}
    scenarios::Vector{Vector{Pair{Int,S}}}
    children::Vector{Vector{Int}} # to children *ba nodes*
    parent_b::Vector{Int} # maps to parent *belief node*
    parent::Vector{Int} # maps to the parent *ba node*
    Delta::Vector{Int}
    mu::Vector{Float64} # needed for ba_mu, excess_uncertainty
    l::Vector{Float64} # needed to select action, excess_uncertainty
    U::Vector{Float64} # needed for blocking
    l_0::Vector{Float64} # needed for find_blocker, backup of l and mu
    obs::Vector{O}

    ba_children::Vector{Vector{Int}}
    ba_mu::Vector{Float64} # needed for next_best
    ba_rho::Vector{Float64} # needed for backup
    ba_Rsum::Vector{Float64} # needed for backup
    ba_action::Vector{A}

    _discount::Float64 # for inferring L in visualization
end

function DESPOT(p::DESPOTPlanner, b_0)
    S = statetype(p.pomdp)
    A = actiontype(p.pomdp)
    O = obstype(p.pomdp)
    root_scenarios = [i=>rand(p.rng, b_0) for i in 1:p.sol.K]

    scenario_belief = ScenarioBelief(root_scenarios, p.rs, 0, b_0)
    L_0, U_0 = bounds(p.bounds, p.pomdp, scenario_belief)

    if p.sol.bounds_warnings
        bounds_sanity_check(p.pomdp, scenario_belief, L_0, U_0)
    end

    return DESPOT{S,A,O}([root_scenarios],
                         [Int[]],
                         [0],
                         [0],
                         [0],
                         [max(L_0, U_0 - p.sol.lambda)],
                         [L_0],
                         [U_0],
                         [L_0],
                         Vector{O}(undef, 1),

                         Vector{Int}[],
                         Float64[],
                         Float64[],
                         Float64[],
                         A[],
                         discount(p.pomdp)
                 )
end

function expand!(D::DESPOT, b::Int, p::DESPOTPlanner)
    S = statetype(p.pomdp)
    A = actiontype(p.pomdp)
    O = obstype(p.pomdp)
    odict = Dict{O, Int}()

    belief = get_belief(D, b, p.rs)
    for a in actions(p.pomdp, belief)
        empty!(odict)
        rsum = 0.0

        for scen in D.scenarios[b]
            rng = get_rng(p.rs, first(scen), D.Delta[b])
            s = last(scen)
            if !isterminal(p.pomdp, s)
                sp, o, r = @gen(:sp, :o, :r)(p.pomdp, s, a, rng)
                rsum += r
                bp = get(odict, o, 0)
                if bp == 0
                    push!(D.scenarios, Vector{Pair{Int, S}}())
                    bp = length(D.scenarios)
                    odict[o] = bp
                end
                push!(D.scenarios[bp], first(scen)=>sp)
            end
        end

        push!(D.ba_children, collect(values(odict)))
        ba = length(D.ba_children)
        push!(D.ba_action, a)
        push!(D.children[b], ba)
        rho = rsum*discount(p.pomdp)^D.Delta[b]/p.sol.K - p.sol.lambda
        push!(D.ba_rho, rho)
        push!(D.ba_Rsum, rsum)

        nbps = length(odict)
        resize!(D, length(D.children) + nbps)
        for (o, bp) in odict
            D.obs[bp] = o
            D.children[bp] = Int[]
            D.parent_b[bp] = b
            D.parent[bp] = ba
            D.Delta[bp] = D.Delta[b]+1

            scenario_belief = get_belief(D, bp, p.rs)
            L_0, U_0 = bounds(p.bounds, p.pomdp, scenario_belief)

            if p.sol.bounds_warnings
                bounds_sanity_check(p.pomdp, scenario_belief, L_0, U_0)
            end

            l_0 = length(D.scenarios[bp])/p.sol.K * discount(p.pomdp)^D.Delta[bp] * L_0
            mu_0 = max(l_0, length(D.scenarios[bp])/p.sol.K * discount(p.pomdp)^D.Delta[bp] * U_0 - p.sol.lambda)

            D.mu[bp] = mu_0
            D.U[bp] = U_0
            D.l[bp] = l_0 # = max(l_0, l_0 - p.sol.lambda)
            D.l_0[bp] = l_0
        end

        push!(D.ba_mu, D.ba_rho[ba] + sum(D.mu[bp] for bp in D.ba_children[ba]))
    end
end

function get_belief(D::DESPOT, b::Int, rs::DESPOTRandomSource)
    if isassigned(D.obs, b)
        ScenarioBelief(D.scenarios[b], rs, D.Delta[b], D.obs[b])
    else
        ScenarioBelief(D.scenarios[b], rs, D.Delta[b], missing)
    end
end

function Base.resize!(D::DESPOT, n::Int)
    resize!(D.children, n)
    resize!(D.parent_b, n)
    resize!(D.parent, n)
    resize!(D.Delta, n)
    resize!(D.mu, n)
    resize!(D.l, n)
    resize!(D.U, n)
    resize!(D.l_0, n)
    resize!(D.obs, n)
end


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

name = "libsndfile"
version = v"1.0.28"

# Collection of sources required to build
sources = [
    ArchiveSource("http://www.mega-nerd.com/libsndfile/files/libsndfile-$(version).tar.gz",
                  "1ff33929f042fa333aed1e8923aa628c3ee9e1eb85512686c55092d1e5a9dfa9")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libsndfile-*/
CFLAGS="-I${prefix}/include" ./configure --prefix=$prefix --host=$target --disable-static
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libsndfile", :libsndfile),
    ExecutableProduct("sndfile-cmp", :sndfile_cmp),
    ExecutableProduct("sndfile-concat", :sndfile_concat),
    ExecutableProduct("sndfile-convert", :sndfile_convert),
    ExecutableProduct("sndfile-deinterleave", :sndfile_deinterleave),
    ExecutableProduct("sndfile-info", :sndfile_info),
    ExecutableProduct("sndfile-interleave", :sndfile_interleave),
    ExecutableProduct("sndfile-metadata-get", :sndfile_metadata_get),
    ExecutableProduct("sndfile-metadata-set", :sndfile_metadata_set),
    ExecutableProduct("sndfile-play", :sndfile_play),
    ExecutableProduct("sndfile-salvage", :sndfile_salvage)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "FLAC_jll",
    "Ogg_jll",
    "libvorbis_jll"
]

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

ext: jl


text:
######################################################################
# General example of how to create & manipulate a window.
# -----
# Licensed under MIT License
using Pkg
Pkg.activate(".")

using KirUtil
using GLFWAbstraction

GLFWAbstraction.on_window_move(::Window{:mywnd}, posx::Integer, posy::Integer) = println("window moved to $posx, $posy")
GLFWAbstraction.on_window_resize(::Window{:mywnd}, sizex::Integer, sizey::Integer) = println("window resized to $sizex×$sizey")
GLFWAbstraction.on_window_close(::Window{:mywnd}) = println("window closing")

@windowhint TransparentFramebuffer true

let wnd = window(:mywnd, "Hello, world!", 960, 480)
    @assert !wnd.shouldclose
    wnd.aspect_ratio = 16//9
    use(wnd)

    while !wnd.shouldclose
        sleep(0.1)
        swapbuffers(wnd)
        pollevents()
    end
end


text:
export Agent

using Base.Threads: @spawn

import Functors

"""
    Agent(;policy, trajectory)

A wrapper of an `AbstractPolicy`. Generally speaking, it does nothing but to
update the trajectory and policy appropriately in different stages.

# Keywords & Fields

- `policy`::[`AbstractPolicy`](@ref): the policy to use
- `trajectory`::[`Trajectory`](@ref): used to store intractions between an agent and an environment
"""
mutable struct Agent{P,T} <: AbstractPolicy
    policy::P
    trajectory::T
    cache::NamedTuple # trajectory do not support partial inserting

    function Agent(policy::P, trajectory::T, cache = NamedTuple()) where {P,T}
        agent = new{P,T}(policy, trajectory, cache)
        if TrajectoryStyle(trajectory) === AsyncTrajectoryStyle()
            bind(trajectory, @spawn(optimise!(p, t)))
        end
        agent
    end
end

Agent(; policy, trajectory, cache = NamedTuple()) = Agent(policy, trajectory, cache)

RLBase.optimise!(agent::Agent) = optimise!(TrajectoryStyle(agent.trajectory), agent)
RLBase.optimise!(::SyncTrajectoryStyle, agent::Agent) =
    optimise!(agent.policy, agent.trajectory)

# already spawn a task to optimise inner policy when initializing the agent
RLBase.optimise!(::AsyncTrajectoryStyle, agent::Agent) = nothing

function RLBase.optimise!(policy::AbstractPolicy, trajectory::Trajectory)
    for batch in trajectory
        optimise!(policy, batch)
    end
end

Functors.functor(x::Agent) =
    (policy = x.policy,), y -> Agent(y.policy, x.trajectory, x.cache)

# !!! TODO: In async scenarios, parameters of the policy may still be updating
# (partially), which will result to incorrect action. This should be addressed
# in Oolong.jl with a wrapper
function (agent::Agent)(env::AbstractEnv)
    action = agent.policy(env)
    push!(agent.trajectory, (agent.cache..., action = action))
    agent.cache = (;)
    action
end

(agent::Agent)(::PreActStage, env::AbstractEnv) =
    agent.cache = (agent.cache..., state = state(env))

(agent::Agent)(::PostActStage, env::AbstractEnv) =
    agent.cache = (agent.cache..., reward = reward(env), terminal = is_terminated(env))


text:
int_rules_5_1_2 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 5.1.2*(d*x)^m*(a+b*arcsin(c*x))^n =#
    #= Antiderivative(((~a')+(~b')*asin((~c')*(~x)))^(~n')/(~x),~x) :=  1/b*Subst(Antiderivative(x^n*cot(-a/b+x/b),x),x,a+b*asin(c*x)) <--  FreeQ([a,b,c],x) && IGtQ(n,0) =#
    #= Antiderivative(((~a')+(~b')*acos((~c')*(~x)))^(~n')/(~x),~x) :=  -1/b*Subst(Antiderivative(x^n*tan(-a/b+x/b),x),x,a+b*acos(c*x)) <--  FreeQ([a,b,c],x) && IGtQ(n,0) =#
    @apply_utils Antiderivative((~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n') / ~x, ~x) => Subst(Antiderivative((~a + ~b * ~x) ^ ~n * cot(~x), ~x), ~x, asin(~c * ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n') / ~x, ~x) => -(Subst(Antiderivative((~a + ~b * ~x) ^ ~n * tan(~x), ~x), ~x, acos(~c * ~x))) <-- FreeQ([~a, ~b, ~c], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => ((~d * ~x) ^ (~m + 1) * (~a + ~b * asin(~c * ~x)) ^ ~n) / (~d * (~m + 1)) - ((~b * ~c * ~n) / (~d * (~m + 1))) * Antiderivative(((~d * ~x) ^ (~m + 1) * (~a + ~b * asin(~c * ~x)) ^ (~n - 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (IGtQ(~n, 0) && NeQ(~m, -1))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => ((~d * ~x) ^ (~m + 1) * (~a + ~b * acos(~c * ~x)) ^ ~n) / (~d * (~m + 1)) + ((~b * ~c * ~n) / (~d * (~m + 1))) * Antiderivative(((~d * ~x) ^ (~m + 1) * (~a + ~b * acos(~c * ~x)) ^ (~n - 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (IGtQ(~n, 0) && NeQ(~m, -1))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~n, ~x) => ((~x) ^ (~m + 1) * (~a + ~b * asin(~c * ~x)) ^ ~n) / (~m + 1) - ((~b * ~c * ~n) / (~m + 1)) * Antiderivative(((~x) ^ (~m + 1) * (~a + ~b * asin(~c * ~x)) ^ (~n - 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (IGtQ(~m, 0) && GtQ(~n, 0))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~n, ~x) => ((~x) ^ (~m + 1) * (~a + ~b * acos(~c * ~x)) ^ ~n) / (~m + 1) + ((~b * ~c * ~n) / (~m + 1)) * Antiderivative(((~x) ^ (~m + 1) * (~a + ~b * acos(~c * ~x)) ^ (~n - 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (IGtQ(~m, 0) && GtQ(~n, 0))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~n, ~x) => ((~x) ^ ~m * sqrt(1 - (~c) ^ 2 * (~x) ^ 2) * (~a + ~b * asin(~c * ~x)) ^ (~n + 1)) / (~b * ~c * (~n + 1)) - (1 / ((~b) ^ 2 * (~c) ^ (~m + 1) * (~n + 1))) * Subst(Antiderivative(ExpandTrigReduce((~x) ^ (~n + 1), sin(-(~a) / ~b + ~x / ~b) ^ (~m - 1) * (~m - (~m + 1) * sin(-(~a) / ~b + ~x / ~b) ^ 2), ~x), ~x), ~x, ~a + ~b * asin(~c * ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IGtQ(~m, 0) && (GeQ(~n, -2) && LtQ(~n, -1)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~n, ~x) => (-((~x) ^ ~m) * sqrt(1 - (~c) ^ 2 * (~x) ^ 2) * (~a + ~b * acos(~c * ~x)) ^ (~n + 1)) / (~b * ~c * (~n + 1)) - (1 / ((~b) ^ 2 * (~c) ^ (~m + 1) * (~n + 1))) * Subst(Antiderivative(ExpandTrigReduce((~x) ^ (~n + 1), cos(-(~a) / ~b + ~x / ~b) ^ (~m - 1) * (~m - (~m + 1) * cos(-(~a) / ~b + ~x / ~b) ^ 2), ~x), ~x), ~x, ~a + ~b * acos(~c * ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (IGtQ(~m, 0) && (GeQ(~n, -2) && LtQ(~n, -1)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~n, ~x) => (((~x) ^ ~m * sqrt(1 - (~c) ^ 2 * (~x) ^ 2) * (~a + ~b * asin(~c * ~x)) ^ (~n + 1)) / (~b * ~c * (~n + 1)) - (~m / (~b * ~c * (~n + 1))) * Antiderivative(((~x) ^ (~m - 1) * (~a + ~b * asin(~c * ~x)) ^ (~n + 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x)) + ((~c * (~m + 1)) / (~b * (~n + 1))) * Antiderivative(((~x) ^ (~m + 1) * (~a + ~b * asin(~c * ~x)) ^ (~n + 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (IGtQ(~m, 0) && LtQ(~n, -2))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~n, ~x) => ((-((~x) ^ ~m) * sqrt(1 - (~c) ^ 2 * (~x) ^ 2) * (~a + ~b * acos(~c * ~x)) ^ (~n + 1)) / (~b * ~c * (~n + 1)) + (~m / (~b * ~c * (~n + 1))) * Antiderivative(((~x) ^ (~m - 1) * (~a + ~b * acos(~c * ~x)) ^ (~n + 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x)) - ((~c * (~m + 1)) / (~b * (~n + 1))) * Antiderivative(((~x) ^ (~m + 1) * (~a + ~b * acos(~c * ~x)) ^ (~n + 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (IGtQ(~m, 0) && LtQ(~n, -2))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~n, ~x) => (1 / (~b * (~c) ^ (~m + 1))) * Subst(Antiderivative((~x) ^ ~n * sin(-(~a) / ~b + ~x / ~b) ^ ~m * cos(-(~a) / ~b + ~x / ~b), ~x), ~x, ~a + ~b * asin(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~n, ~x) => (-1 / (~b * (~c) ^ (~m + 1))) * Subst(Antiderivative((~x) ^ ~n * cos(-(~a) / ~b + ~x / ~b) ^ ~m * sin(-(~a) / ~b + ~x / ~b), ~x), ~x, ~a + ~b * acos(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~n], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => Unintegrable((~d * ~x) ^ ~m * (~a + ~b * asin(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * (~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => Unintegrable((~d * ~x) ^ ~m * (~a + ~b * acos(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x)
end

max_stars_repo_path: src/rules/5 Inverse trig functions/5.1 Inverse sine/5.1.2 (d x)^m (a+b arcsin(c x))^n.jl


text:
const RemotePool = PoolType{Remote}
const LocalPool = PoolType{Local}

#####
##### Backedges
#####

# Back edges from the manager to the GC managed object holding a block.
# Kept as a raw Ptr{Ptr} to avoid GC cycles.
const Backedge = Ptr{Ptr{Nothing}}
unsafe_block(backedge::Backedge) = unsafe_block(unsafe_load(backedge))

backedge(x::Backedge) = x
backedge(x::Ptr) = convert(Backedge, x)
backedge(x) = backedge(blockpointer(x))

mutable struct BackedgeMap
    dict::Dict{UInt,Tuple{Block,Backedge}}
    size::Int
end

BackedgeMap() = BackedgeMap(Dict{UInt,Backedge}(), 0)

Base.getindex(map::BackedgeMap, id) = map.dict[id]
Base.setindex!(map::BackedgeMap, bb::Tuple{Block,Backedge}, id) = setindex!(map.dict, bb, id)
function set!(map::BackedgeMap, block::Block, backedge::Backedge)
    id = getid(block)
    sz = length(block)

    @check !in(id, map)
    map.dict[id] = (block, backedge)
    map.size += sz
    return backedge
end

function Base.delete!(map::BackedgeMap, id::UInt, sz::Integer)
    delete!(map.dict, id)
    @check map.size >= sz
    map.size -= sz
    return nothing
end

getsize(map::BackedgeMap) = map.size
Base.in(id, map::BackedgeMap) = haskey(map.dict, id)
Base.haskey(map::BackedgeMap, id) = haskey(map.dict, id)
Base.length(map::BackedgeMap) = length(map.dict)
Base.isempty(map::BackedgeMap) = isempty(map.dict)
Base.keys(map::BackedgeMap) = keys(map.dict)
Base.values(map::BackedgeMap) = values(map.dict)

Base.iterate(map::BackedgeMap, s...) = iterate(map.dict, s...)

#####
##### Cache Manager
#####

struct NoTelemetry end
mutable struct CacheManager{C,T}
    map::BackedgeMap

    # local datastructures
    policy::C
    remote_heap::CompactHeap
    local_heap::CompactHeap

    # Create a new ID for each object registered in the cache.
    idcount::Threads.Atomic{UInt64}

    # Synchronize access
    alloc_lock::Base.Threads.SpinLock
    freebuffer::FreeBuffer{Block}
    telemetry::T
end

# Telemetry Utils
gettelemetry(manager::CacheManager) = manager.telemetry
telemetry_enabled(::CacheManager) = true
telemetry_enabled(::CacheManager{<:Any,NoTelemetry}) = false
getpolicy(manager::CacheManager) = manager.policy

macro telemetry(manager, expr)
    return quote
        if !isa($(esc(manager)).telemetry, NoTelemetry)
            $(esc(expr))
        end
    end
end

# Accessor Utils
getmap(manager::CacheManager, id::Integer) = manager.map[id]
getmap(manager::CacheManager) = manager.map

getheap(manager::CacheManager, ::PoolType{Local}) = manager.local_heap
getheap(manager::CacheManager, ::PoolType{Remote}) = manager.remote_heap

candrain(manager::CacheManager) = candrain(manager.freebuffer)

"""
$(TYPEDSIGNATURES)

Returns the allocation lock for `manager` without acquiring it.
"""
alloc_lock(manager::CacheManager) = manager.alloc_lock

mb(x) = x * 1_000_000

function CacheManager(
    local_allocator,
    remote_allocator;
    localsize = mb(1_000),
    remotesize = mb(1_000),
    minallocation = 10,
    policy = OptaneTracker((2^minallocation,)),
    telemetry = NoTelemetry(),
)
    map = BackedgeMap()

    # Initialize Heaps
    remote_heap = CompactHeap(remote_allocator, remotesize; pool = Remote, minallocation)
    local_heap = CompactHeap(local_allocator, localsize; pool = Local, minallocation)

    # Construct the manager.
    manager = CacheManager(
        map,
        policy,
        remote_heap,
        local_heap,
        Threads.Atomic{UInt64}(1),
        Threads.SpinLock(),
        FreeBuffer{Block}(),
        telemetry,
    )

    # Add this to the global manager list to ensure that it outlives any of its users.
    push!(GlobalManagers, manager)
    return manager
end

getid(manager::CacheManager) = Threads.atomic_add!(manager.idcount, one(UInt64))
readid(manager::CacheManager) = manager.idcount[]

# Can only GC if all objects tracked via backedges have been collected.
function cangc(manager::CacheManager)
    # Empty out the cleanlist to deal with anything that's been GC'd
    @spinlock alloc_lock(manager) unsafe_cleanup!(manager)
    return isempty(manager.map)
end

function Base.show(io::IO, M::CacheManager)
    println(io, "Cache Manager")
    println(io, "    $(length(M.map)) Objects")
    println(io, "    $(getsize(M.map) / 1E9) GB Memory Used.")
    allocated, total = getstate(M.local_heap) ./ 1E9
    println(io, "    Local Heap utilization: $allocated of $total ($(allocated / total) %)")
    allocated, total = getstate(M.remote_heap) ./ 1E9
    println(
        io,
        "    Remote Heap utilization: $allocated of $total ($(allocated / total) %)",
    )
end

#####
##### Object (De)Registration
#####

function unsafe_register!(manager::CacheManager, object::Object)
    @requires alloc_lock(manager)
    block = metadata(object)
    set!(getmap(manager), block, backedge(object))
    push!(manager.policy, block)
    return nothing
end

function unsafe_unregister!(manager::CacheManager, block::Block)
    @requires alloc_lock(manager)
    delete!(getmap(manager), getid(block), length(block))
    delete!(manager.policy, block)
    return nothing
end

#####
##### Lowest Level Alloc and Free
#####

function unsafe_free_direct(manager::CacheManager, block::Block)
    ptr = datapointer(block)
    pool = getpool(block)
    @telemetry manager telemetry_dealloc(gettelemetry(manager), block.id, pool)
    if pool == Local
        free(getheap(manager, LocalPool()), ptr)
    elseif pool == Remote
        free(getheap(manager, RemotePool()), ptr)
    end
    return nothing
end

function unsafe_alloc_direct(
    pool::PoolType{T}, manager::CacheManager, bytes, id::UInt
) where {T}
    @requires alloc_lock(manager)
    ptr = alloc(getheap(manager, pool), bytes, id)
    if ptr !== nothing
        @telemetry manager telemetry_alloc(gettelemetry(manager), bytes, id, ptr, T)
    end
    return ptr
end

"""
    direct_alloc(pool, manager, bytes, id) -> Object

Allocate an object from `manager` in `pool`.
"""
function direct_alloc(pool::PoolType, manager::CacheManager, bytes, id::UInt)
    @spinlock alloc_lock(manager) begin
        ptr = iszero(bytes) ? Ptr{Nothing}() : unsafe_alloc_direct(pool, manager, bytes, id)
        ptr === nothing && throw(AllocationException())
        return Object(ptr, manager)
    end
end

#####
##### Cleanup Finalizer
#####

# In general, we don't know when the GC will run.
# So, we make the finalizer really short, just appending the block to be cleaned up
# to a vector.
#
# When we're trying to allocate (i.e., holding the lock) - THEN we'll call the `_cleanup`
# method below which will put back all of the blocks on the `cleanlist`.
prepare_cleanup!(manager::CacheManager) = unsafe_swap!(manager.freebuffer)
maybe_cleanup!(M::CacheManager, args...) = candrain(M) && unsafe_cleanup!(M, args...)
function unsafe_cleanup!(M::CacheManager, id = nothing)
    @requires alloc_lock(M)
    prepare_cleanup!(M)
    id_cleaned = false

    while true
        cleanlist = unsafe_get(M.freebuffer)

        # Free all blocks in the cleanlist
        for block in cleanlist
            @check !isfree(block)

            id == getid(block) && (id_cleaned = true)
            @telemetry M telemetry_gc(gettelemetry(M), getid(block))
            sibling = getsibling(block)
            if sibling !== nothing
                @check getid(sibling) == getid(block)
                @check getpool(sibling) != getpool(block)
                @check getsibling(sibling) == block

                delete!(M.policy, sibling)
                unsafe_free_direct(M, sibling)
            end
            unsafe_unregister!(M, block)
            unsafe_free_direct(M, block)
        end
        empty!(cleanlist)

        # Process other buffer as well.
        # Safe to do since we already hold the "remove_lock".
        # The next trip around the loop will be with the new buffer.
        if candrain(M.freebuffer)
            unsafe_swap!(M.freebuffer)
        else
            break
        end
    end
    return id_cleaned
end

#####
##### Policy Level Allocations
#####

defrag!(manager::CacheManager) = defrag!(manager, manager.policy)

function unsafe_alloc_through_policy(
    manager::CacheManager,
    bytes,
    priority::AllocationPriority,
    id = getid(manager),
)
    # Only attempt draining at the beginning of an allocation round.
    #safeprint("Allocating $(bytes) for id $(id)"; force = true)
    candrain(manager.freebuffer) && unsafe_cleanup!(manager)
    ptr = policy_new_alloc(manager.policy, manager, bytes, id, priority)
    ptr === nothing && error("Allocation failed!")
    return ptr
end

# Alloc through the policy
function alloc(
    manager::CacheManager,
    bytes::Int,
    priority::AllocationPriority = PreferLocal,
    id::UInt = getid(manager),
)
    @spinlock alloc_lock(manager) begin
        # TODO: What to do about zero sized allocations ...
        ptr =
            iszero(bytes) ? Ptr{Nothing}() :
            unsafe_alloc_through_policy(manager, bytes, priority, id)
        return Object(ptr, manager)
    end
end

#####
##### Top level Frees
#####

# Note: `push!` for the `FreeBuffer` is thread safe, so no synchonizations needs
# to happen at this level.
# @static if ALLOW_UNSAFE_FREE

# If we allow "unsafe_free", then the pointer in a object will be replaced by
# a null pointer, in which case we definitely don't want to add this to the
# free buffer.
function free(manager::CacheManager, block::Block)
    push!(manager.freebuffer, block)
end
function free(manager::CacheManager, ptr::Ptr)
    isnull(ptr) || free(manager, unsafe_block(ptr))
end

function unsafe_free(object::Object)
    ptrptr = Ptr{Ptr{Nothing}}(blockpointer(object))
    old = atomic_ptr_xchg!(ptrptr, Ptr{Nothing}())
    if !isnull(old)
        block = unsafe_block(old)
        manager = object.manager
        @telemetry manager telemetry_gc(gettelemetry(manager), getid(block))
        free(manager, block)
    end
end

function unsafe_free(manager::CacheManager, (block, ptrptr)::Tuple{Block,Backedge})
    isqueued(block) && return nothing
    old = atomic_ptr_xchg!(ptrptr, Ptr{Nothing}())
    if !isnull(old)
        @telemetry manager telemetry_unsafefree(gettelemetry(manager), getid(block))
        free(manager, block)
    end
    return nothing
end

# else
#     free(manager::CacheManager, ptr::Ptr) = push!(manager.freebuffer, unsafe_block(ptr))
#     unsafe_free(::Object) = nothing
# end

#####
##### Set primary
#####

"""
    unsafe_setprimary!(manager::CacheManager, current::Block, next::Block)::Union{Block, Nothing}

Switch a primary segment from `current` to `next`.
This operation may fail if `current` has been garbage collected.
If this is the case, this function will return `nothing`.
"""
function unsafe_setprimary!(
    manager::CacheManager,
    current::Block,
    next::Block;
    unsafe = false,
)
    # Need to check if the current is queued for cleanup.
    # If so, this means that the backedge is invalid.
    isqueued(current) && return nothing
    id = getid(current)
    unsafe || @check getid(next) == id
    map = getmap(manager)
    primary, backedge = map[id]
    @check primary == current

    # Small chance GC ran. Make sure that `current` is still not queued before comitting.
    isqueued(primary) && return nothing
    old = atomic_ptr_xchg!(backedge, datapointer(next))
    map[id] = (next, backedge)
    unsafe || @telemetry manager telemetry_primary(gettelemetry(manager), id, getpool(next))

    # Check that in fact the passed `current` block IS the primary segment.
    @check old == datapointer(current)
    return old
end

function unsafe_clearprimary!(manager::CacheManager, current::Block)
    isqueued(current) && return nothing
    id = getid(current)
    map = getmap(manager)

    _, backedge = map[id]
    old = atomic_ptr_xchg!(backedge, Ptr{Nothing}())
    map[id] = (Block(), backedge)

    @check old == datapointer(current)
    return old
end

getprimary(manager::CacheManager, block::Block) = getprimary(manager, getid(block))
getprimary(manager::CacheManager, id::UInt) = first(getmap(manager, id))

#####
##### Copyto!
#####

function Base.copyto!(dst::Block, src::Block, ::CacheManager; include_header = false)
    if src.size <= 4096
        nthreads = 1
    else
        nthreads = getpool(dst) == Remote ? 8 : Threads.nthreads()
    end

    if include_header
        _memcpy!(pointer(dst), pointer(src), length(src); nthreads)
    else
        _memcpy!(datapointer(dst), datapointer(src), length(src); nthreads)
    end
    return nothing
end

#####
##### Utils
#####

function visible_ids(manager::CacheManager, pool; primary_only = false)
    ids = Set{UInt64}()
    for (id, bb) in getmap(manager)
        _, ptr = bb
        block = unsafe_block(ptr)
        sibling = getsibling(block)
        if getpool(block) == pool
            push!(ids, id)
        elseif !primary_only && (sibling !== nothing)
            @check getpool(sibling) == pool
            push!(ids, id)
        end
    end
    return ids
end

#####
##### Policy
#####

@inline function setdirty!(manager::CacheManager, block::Block, flag = true)
    return setdirty!(manager.policy, block, flag)
end

update!(manager::CacheManager, block::Block) = update!(manager.policy, block)


text:
module FastSpike

using LinearAlgebra
using Random
using RecursiveArrayTools
using JLD2
using NNlib
using CUDA
using Adapt
using OMEinsum
using Parameters: @with_kw

export add_group!, connect!, run!, train!, reset!, pad1D, pad2D, record!, get, save, load,
        gpu, cpu, randomConnection, regularActivity, randomActivity, makeInput, EI,
        raster, set_reward
export SpikingNetwork, Network, DelayNetwork, LearningRule, NeuronType, NeuronGroup, LIF,
        Izhikevich, STDP, Monitor, Connection, UpdateRule, RewardModulatedUpdateRule,
        RegularUpdate, WeightDependent, WeightDependentUpdate, Softbound, RewardModulated,
        WeightDependentRewardModulated

abstract type SpikingNetwork{T} end

include("utils.jl")
include("connection.jl")
include("neuron.jl")
include("update_rules.jl")
include("learning_rules.jl")
include("network.jl")
include("delay_network.jl")
include("monitor.jl")
include("learning.jl")
include("cuda.jl")

end


text:
function calcx(E, m, n, k)
    tmp = E.'*E;
    R = E * ( tmp - spdiagm( diag(tmp) ) );
    r,c,v = findnz(R);
    id = v.==2;
    A = sparse( r[id], c[id], 1, m, n);
    s = sum(A, 2);
    x = s .< (k-2);
    return(x, !x);
end

function ktruss(inc_mtx_file, k)
    if ~isfile( inc_mtx_file )
        println("unable to open input file");
        return (-1);
    end

    # load input data
    t_read_inc=@elapsed ii = readdlm( inc_mtx_file, '\t', Int64);
    println("incidence matrix read time : ", t_read_inc);

    t_create_inc=@elapsed E = sparse( ii[:,1], ii[:,2], ii[:,3] );
    println("sparse adj. matrix creation time : ", t_create_inc);

    #
    tic();
    m,n = size(E);
    x, xc = calcx(E, m, n, k);
    while sum(xc) != sum( any(E,2) )
        E[find(x), :] = 0
        x, xc = calcx(E, m, n, k);
    end

    return E
end

#######################################################
# Graph Challenge Benchmark
# Architect : Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Developer : Dr. Siddharth Samsi (ssamsi@mit.edu)
#
# MIT
########################################################
# (c) <2017> Massachusetts Institute of Technology
########################################################

max_stars_repo_name: bkmgit/GraphChallenge

max_issues_repo_name: bkmgit/GraphChallenge

max_forks_repo_name: bkmgit/GraphChallenge


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

export
    sssr, atom_sssrmem, bond_sssrmem,
    atom_sssrsizes, bond_sssrsizes,
    isringatom, isringbond,
    atom_sssrcount, bond_sssrcount,
    scaffoldmem, componentmem,
    atomsymbol, charge, multiplicity, bondorder, nodedegree,
    apparentvalence, valence, lonepair, heavyatomcount,
    explicithcount, implicithcount, hcount,
    connectivity, pielectron, ishdonor, ishacceptor, stdweight,
    isrotatable, isaromatic, isaromaticbond,
    molweight,
    hacceptorcount,
    hdonorcount,
    wclogp,
    rotatablecount

# Molecular graph topology

sssr = mincyclenodes
atom_sssrmem = node_cyclemem
bond_sssrmem = edge_cyclemem
scaffoldmem = two_edge_membership
componentmem = connected_membership

function atom_sssrsizes(mol::UndirectedGraph)
    vec = [Set{Int}() for i in 1:nodecount(mol)]
    rings = sssr(mol)
    for (i, cs) in enumerate(atom_sssrmem(mol))
        for c in cs
            push!(vec[i], length(rings[c]))
        end
    end
    return vec
end

atom_sssrsizes(view::SubgraphView) = atom_sssrsizes(view.graph)

function bond_sssrsizes(mol::UndirectedGraph)
    vec = [Set{Int}() for i in 1:edgecount(mol)]
    rings = sssr(mol)
    for (i, cs) in enumerate(bond_sssrmem(mol))
        for c in cs
            push!(vec[i], length(rings[c]))
        end
    end
    return vec
end

bond_sssrsizes(view::SubgraphView) = bond_sssrsizes(view.graph)

# TODO: waiting for fix #28992
# @cache atom_isringmem(mol) = .!isempty.(node_cyclemem(mol))
isringatom(mol::UndirectedGraph) = [!i for i in isempty.(atom_sssrmem(mol))]
isringatom(view::SubgraphView) = isringatom(view.graph)

# @cache bond_isringmem(mol) = .!isempty.(node_cyclemem(mol))
isringbond(mol::UndirectedGraph) = [!i for i in isempty.(bond_sssrmem(mol))]
isringbond(view::SubgraphView) = isringbond(view.graph)

atom_sssrcount(mol::UndirectedGraph) = length.(atom_sssrmem(mol))
atom_sssrcount(view::SubgraphView) = atom_sssrcount(view.graph)

bond_sssrcount(mol::UndirectedGraph) = length.(bond_sssrmem(mol))
bond_sssrcount(view::SubgraphView) = bond_sssrcount(view.graph)

# Elemental properties

"""
    atomsymbol(mol::GraphMol) -> Vector{Int}

Return the size `n` vector of atom symbols within the molecule that have `n`
atoms.
"""
@cache atomsymbol(mol::GraphMol) = getproperty.(nodeattrs(mol), :symbol)
atomsymbol(view::SubgraphView) = atomsymbol(view.graph)

"""
    charge(mol::GraphMol) -> Vector{Int}

Return the size `n` vector of atom charges within the molecule that have `n`
atoms.
"""
@cache charge(mol::GraphMol) = getproperty.(nodeattrs(mol), :charge)
charge(view::SubgraphView) = charge(view.graph)

"""
    multiplicity(mol::GraphMol) -> Vector{Int}

Return the size `n` vector of the atom multiplicity within the molecule that
have `n` atoms (1: non-radical, 2: radical, 3: biradical).
"""
@cache multiplicity(mol::GraphMol) = getproperty.(nodeattrs(mol), :multiplicity)
multiplicity(view::SubgraphView) = multiplicity(view.graph)

"""
    bondorder(mol::GraphMol) -> Vector{Int}

Return the size `n` vector of the bond order within the molecule that have `n`
bonds.
"""
@cache bondorder(mol::GraphMol) = getproperty.(edgeattrs(mol), :order)
bondorder(view::SubgraphView) = bondorder(view.graph)

"""
    nodedegree(mol::GraphMol) -> Vector{Int}

Return the size `n` vector of the node degree within the molecular graph that
have `n` atom nodes.
"""
@cache nodedegree(mol::GraphMol) = [degree(mol, n) for n in 1:nodecount(mol)]
nodedegree(view::SubgraphView) = nodedegree(view.graph)

"""
    apparentvalence(mol::GraphMol) -> Vector{Int}

Return the size `n` vector of the apparent valence (the sum of incident bond
order) within the molecule that have `n` atoms.
"""
@cache function apparentvalence(mol::GraphMol)
    vec = zeros(Int, nodecount(mol))
    bondorder_ = bondorder(mol)
    for i in 1:nodecount(mol)
        for inc in incidences(mol, i)
            vec[i] += bondorder_[inc]
        end
    end
    return vec
end

apparentvalence(view::SubgraphView) = apparentvalence(view.graph)

"""
    valence(mol::GraphMol) -> Vector{Union{Int,Nothing}}

Return the size `n` vector of the intrinsic valence (with considering implicit
hydrogens) within the molecule that have `n` atoms.

Note that implicit hydrogens are available for only organic atoms. The intrinsic
valence value of inorganic atoms would be `nothing`.
"""
@cache function valence(mol::GraphMol)
    defs = Dict(
        :H => 1, :B => 3, :C => 4, :N => 3, :O => 2, :F => 1,
        :Si => 4, :P => 3, :S => 2, :Cl => 1,
        :As => 3, :Se => 2, :Br => 1, :I => 1
    )
    vec = Union{Int,Nothing}[]
    for (sym, chg) in zip(atomsymbol(mol), charge(mol))
        num = get(defs, sym, nothing)
        v = num === nothing ? nothing : num + chg
        push!(vec, v)
    end
    return vec
end

valence(view::SubgraphView) = valence(view.graph)

"""
    lonepair(mol::GraphMol) -> Vector{Union{Int,Nothing}}

Return the size `n` vector of the number of lone pairs within the molecule that
have `n` atoms.

Note that implicit hydrogens are available for only organic atoms. The lonepair
value of inorganic atoms would be `nothing`.
"""
@cache function lonepair(mol::GraphMol)
    defs = Dict(
        :H => 0, :B => -1, :C => 0, :N => 1, :O => 2, :F => 3,
        :Si => 0, :P => 1, :S => 2, :Cl => 3,
        :As => 1, :Se => 2, :Br => 3, :I => 3
    )
    vec = Union{Int,Nothing}[]
    for (sym, chg) in zip(atomsymbol(mol), charge(mol))
        num = get(defs, sym, nothing)
        v = num === nothing ? nothing : num - chg
        push!(vec, v)
    end
    return vec
end

lonepair(view::SubgraphView) = lonepair(view.graph)

"""
    heavyatomcount(mol::GraphMol) -> Vector{Int}

Return the size `n` vector of the number of adjacent non-hydrogen atoms within
the molecule that have `n` atoms.
"""
@cache function heavyatomcount(mol::GraphMol)
    vec = zeros(Int, nodecount(mol))
    atomsymbol_ = atomsymbol(mol)
    for i in 1:nodecount(mol)
        for adj in adjacencies(mol, i)
            if atomsymbol_[adj] != :H
                vec[i] += 1
            end
        end
    end
    return vec
end

heavyatomcount(view::SubgraphView) = heavyatomcount(view.graph)

@cache function implicithcount(mol::GraphMol)
    hcnt = (v, av) -> v === nothing ? 0 : max(0, v - av)
    return hcnt.(valence(mol), apparentvalence(mol))
end

implicithcount(view::SubgraphView) = implicithcount(view.graph)

@cache explicithcount(mol::GraphMol) = nodedegree(mol) - heavyatomcount(mol)
explicithcount(view::SubgraphView) = explicithcount(view.graph)

@cache hcount(mol::GraphMol) = explicithcount(mol) + implicithcount(mol)
hcount(view::SubgraphView) = hcount(view.graph)

@cache connectivity(mol::GraphMol) = nodedegree(mol) + implicithcount(mol)
connectivity(view::SubgraphView) = connectivity(view.graph)

@cache pielectron(mol::GraphMol) = apparentvalence(mol) - nodedegree(mol)
pielectron(view::SubgraphView) = pielectron(view.graph)

@cache function ishdonor(mol::GraphMol)
    dc = (sym, h) -> sym in (:N, :O) && h > 0
    return dc.(atomsymbol(mol), hcount(mol))
end

ishdonor(view::SubgraphView) = ishdonor(view.graph)

@cache function ishacceptor(mol::GraphMol)
    ac = (sym, lp) -> lp === nothing ? false : sym in (:N, :O, :F) && lp > 0
    return ac.(atomsymbol(mol), lonepair(mol))
end

ishacceptor(view::SubgraphView) = ishacceptor(view.graph)

function stdweight(mol::GraphMol)
    weight = (atom, imh) -> atomweight(atom) + H_WEIGHT * imh
    return weight.(nodeattrs(mol), implicithcount(mol))
end

stdweight(view::SubgraphView) = stdweight(view.graph)

# Rotatable bonds

"""
    isrotatable(mol::GraphMol)

Return whether the bonds are rotatable or not.
"""
@cache function isrotatable(mol::GraphMol)
    nodedegree_ = nodedegree(mol)
    isringbond_ = isringbond(mol)
    bondorder_ = bondorder(mol)
    vec = Bool[]
    for (i, (u, v)) in enumerate(edgesiter(mol))
        rot = (!isringbond_[i] && bondorder_[i] == 1
            && nodedegree_[u] != 1 && nodedegree_[v] != 1)
        push!(vec, rot)
    end
    return vec
end

isrotatable(view::SubgraphView) = isrotatable(view.graph)

# Aromatic rings

@cache function isaromaticring(mol::GraphMol)
    vec = falses(circuitrank(mol))
    for (i, ring) in enumerate(sssr(mol))
        if satisfyhuckel(mol, ring)
            vec[i] = true
        elseif nodeattrtype(mol) === SmilesAtom # SMILES aromatic atom
            sub = nodesubgraph(mol, Set(ring))
            if all(nodeattr(mol, n).isaromatic for n in nodeset(sub))
                vec[i] = true
            end
        end
    end
    return vec
end

"""
    isaromatic(mol::GraphMol)

Return the vector whether the atom belongs to an aromatic ring.

Note that aromaticity described here means simplified binary descriptor
(aromatic or not) based on classical Huckel's rule. This is intended for use in
some kind of pharmaceutical researches. Non-classical aromaticity such as
Moebius aromaticity is not considered.
"""
@cache function isaromatic(mol::GraphMol)
    aromatic = falses(nodecount(mol))
    for ring in sssr(mol)[isaromaticring(mol)]
        sub = nodesubgraph(mol, Set(ring))
        for n in nodeset(sub)
            aromatic[n] = true
        end
    end
    return aromatic
end

isaromatic(view::SubgraphView) = isaromatic(view.graph)

@cache function isaromaticbond(mol::GraphMol)
    aromaticbond = falses(edgecount(mol))
    for ring in sssr(mol)[isaromaticring(mol)]
        sub = nodesubgraph(mol, Set(ring))
        for e in edgeset(sub)
            aromaticbond[e] = true
        end
    end
    return aromaticbond
end

isaromaticbond(view::SubgraphView) = isaromaticbond(view.graph)

function satisfyhuckel(mol::GraphMol, ring)
    cnt = 0
    atomsymbol_ = atomsymbol(mol)
    nodedegree_ = nodedegree(mol)
    pielectron_ = pielectron(mol)
    lonepair_ = lonepair(mol)
    carbonylO = findall(
        (atomsymbol_ .== :O) .* (nodedegree_ .== 1) .* (pielectron_ .== 1))
    carbonylC = Int[]
    for o in carbonylO
        c = iterate(adjacencies(mol, o))[1]
        if atomsymbol_[c] == :C
            push!(carbonylC, c)
        end
    end
    for r in ring
        if r in carbonylC
            continue
        elseif pielectron_[r] == 1
            cnt += 1
        elseif lonepair_[r] === nothing
            return false
        elseif lonepair_[r] > 0
            cnt += 2
        elseif lonepair_[r] < 0
            continue
        else
            return false
        end
    end
    return cnt % 4 == 2
end

# Molecular properties

"""
    molweight(mol::GraphMol; digits=2) -> Float64

Return standard molecular weight.
"""
molweight(mol::GraphMol; digits=2
    ) = round(reduce(+, stdweight(mol); init=0), digits=digits)

"""
    hacceptorcount(mol::GraphMol) -> Int

Return the number of hydrogen bond acceptors (N, O and F).
"""
hacceptorcount(mol::GraphMol) = reduce(+, ishacceptor(mol); init=0)

"""
    hdonorcount(mol::GraphMol) -> Int

Return the number of hydrogen bond donors (O and N attached to hydrogens).
"""
hdonorcount(mol::GraphMol) = reduce(+, ishdonor(mol); init=0)

"""
    wclogp(mol::GraphMol) -> Float64

Return predicted logP value calculated by using Wildman and Crippen method.

# Reference

1. Wildman, S. A. and Crippen, G. M. (1999). Prediction of Physicochemical
   Parameters by Atomic Contributions. Journal of Chemical Information and
   Modeling, 39(5), 868–873. https://doi.org/10.1021/ci990307l
"""
wclogp(mol::GraphMol; digits=2
    ) = round(reduce(+, wclogpcontrib(mol); init=0), digits=digits)

"""
    rotatablecount(mol::GraphMol) -> Int

Return the number of rotatable bonds.
"""
rotatablecount(mol::GraphMol) = reduce(+, isrotatable(mol); init=0)


text:
function get_followers_ids(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/followers/ids.json", options))
end

function get_followers_list(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/followers/list.json", options))
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
"""
$(SIGNATURES)

Calculate slope of particle distribution function `a = a1 + u * a2 + 0.5 * u^2 * a3`
"""
pdf_slope(u, Δ) = Δ / (u + 1e-7)

"""
$(SIGNATURES)
"""
function pdf_slope(prim::AV, sw::AV, inK)
    sl = similar(sw, axes(prim))

    if length(prim) == 3

        sl[3] =
            4.0 * prim[3]^2 / (inK + 1.0) / prim[1] * (
                2.0 * sw[3] - 2.0 * prim[2] * sw[2] +
                sw[1] * (prim[2]^2 - 0.5 * (inK + 1.0) / prim[3])
            )
        sl[2] = 2.0 * prim[3] / prim[1] * (sw[2] - prim[2] * sw[1]) - prim[2] * sl[3]
        sl[1] =
            sw[1] / prim[1] - prim[2] * sl[2] -
            0.5 * (prim[2]^2 + 0.5 * (inK + 1.0) / prim[3]) * sl[3]

    elseif length(prim) == 4

        sl[4] =
            4.0 * prim[4]^2 / (inK + 2.0) / prim[1] * (
                2.0 * sw[4] - 2.0 * prim[2] * sw[2] - 2.0 * prim[3] * sw[3] +
                sw[1] * (prim[2]^2 + prim[3]^2 - 0.5 * (inK + 2.0) / prim[4])
            )
        sl[3] = 2.0 * prim[4] / prim[1] * (sw[3] - prim[3] * sw[1]) - prim[3] * sl[4]
        sl[2] = 2.0 * prim[4] / prim[1] * (sw[2] - prim[2] * sw[1]) - prim[2] * sl[4]
        sl[1] =
            sw[1] / prim[1] - prim[2] * sl[2] - prim[3] * sl[3] -
            0.5 * (prim[2]^2 + prim[3]^2 + 0.5 * (inK + 2.0) / prim[4]) * sl[4]

    elseif length(prim) == 5

        sl[5] =
            4.0 * prim[5]^2 / (inK + 3.0) / prim[1] * (
                2.0 * sw[5] - 2.0 * prim[2] * sw[2] - 2.0 * prim[3] * sw[3] -
                2.0 * prim[4] * sw[4] +
                sw[1] * (prim[2]^2 + prim[3]^2 + prim[4]^2 - 0.5 * (inK + 3.0) / prim[5])
            )
        sl[4] = 2.0 * prim[5] / prim[1] * (sw[4] - prim[4] * sw[1]) - prim[4] * sl[5]
        sl[3] = 2.0 * prim[5] / prim[1] * (sw[3] - prim[3] * sw[1]) - prim[3] * sl[5]
        sl[2] = 2.0 * prim[5] / prim[1] * (sw[2] - prim[2] * sw[1]) - prim[2] * sl[5]
        sl[1] =
            sw[1] / prim[1] - prim[2] * sl[2] - prim[3] * sl[3] - prim[4] * sl[4] -
            0.5 * (prim[2]^2 + prim[3]^2 + prim[4]^2 + 0.5 * (inK + 3.0) / prim[5]) * sl[5]

    end

    return sl
end

"""
$(SIGNATURES)

Calculate slope of multi-component particle distribution function `a = a1 + u * a2 + 0.5 * u^2 * a3`

"""
function mixture_pdf_slope(prim::AM, sw::AM, inK)
    sl = similar(sw, axes(prim))

    for j in axes(sl, 2)
        sl[:, j] .= pdf_slope(prim[:, j], sw[:, j], inK)
    end

    return sl
end

"""
$(SIGNATURES)

Reduced distribution function

* @arg : particle distribution function with full velocity space
* @arg : quadrature weights with reduced velocity setting (v & w by default)
"""
function reduce_distribution(f::AM, weights::AV, dim = 1)
    if dim == 1
        h = similar(f, axes(f, 1))
        for i in eachindex(h)
            h[i] = sum(@. weights * f[i, :])
        end
    elseif dim == 2
        h = similar(f, axes(f, 2))
        for j in eachindex(h)
            h[j] = sum(@. weights * f[:, j])
        end
    else
        throw("dimension dismatch")
    end

    return h
end

"""
$(SIGNATURES)
"""
function reduce_distribution(f::AA{T1,3}, weights::AM, dim = 1) where {T1}
    if dim == 1
        h = similar(f, axes(f, 1))
        for i in eachindex(h)
            h[i] = sum(@. weights * f[i, :, :])
        end
    elseif dim == 2
        h = similar(f, axes(f, 2))
        for j in eachindex(h)
            h[j] = sum(@. weights * f[:, j, :])
        end
    elseif dim == 3
        h = similar(f, axes(f, 3))
        for k in eachindex(h)
            h[k] = sum(@. weights * f[:, :, k])
        end
    else
    end

    return h
end

"""
$(SIGNATURES)
"""
function reduce_distribution(
    f::AA{X,3},
    v::Y,
    w::Y,
    weights::AM,
    dim = 1,
) where {X,Y<:AA{<:Real,3}}

    if dim == 1
        h = similar(f, axes(f, 1))
        b = similar(h)
        for i in eachindex(h)
            h[i] = sum(@. weights * f[i, :, :])
            b[i] = sum(@. weights * (v[i, :, :]^2 + w[i, :, :]^2) * f[i, :, :])
        end
    elseif dim == 2
        h = similar(f, axes(f, 2))
        b = similar(h)
        for j in eachindex(h)
            h[j] = sum(@. weights * f[:, j, :])
            b[j] = sum(@. weights * (v[:, j, :]^2 + w[:, j, :]^2) * f[:, j, :])
        end
    elseif dim == 3
        h = similar(f, axes(f, 3))
        b = similar(h)
        for k in eachindex(h)
            h[k] = sum(@. weights * f[:, :, k])
            b[k] = sum(@. weights * (v[:, :, k]^2 + w[:, :, k]^2) * f[:, :, k])
        end
    end

    return h, b

end

"""
$(SIGNATURES)

Recover full distribution function from reduced ones

* @arg h & b : reduced particle distribution function with 1D velocity space
* @arg u : quadrature nodes in 1D velocity space
* @arg weights : quadrature weights in 1D velocity space
* @arg v & w : quadrature nodes in the rest velocity space (with 3D setting)
* @return f : particle distribution function with 3D velocity space
"""
function full_distribution(
    h::X,
    b::X,
    u::Y,
    weights::Y,
    v::Z,
    w::Z,
    ρ,
    γ = 5 / 3,
) where {X<:AA{<:FN,1},Y<:AA{<:FN,1},Z<:AA{<:FN,3}}

    @assert length(h) == size(v, 1) throw(
        DimensionMismatch("reduced and full distribution function mismatch"),
    )

    Ei = 0.5 * discrete_moments(b, u, weights, 0)
    λi = 0.5 * ρ / (γ - 1.0) / Ei / 3.0 * 2.0

    f = similar(v)
    for k in axes(f, 3), j in axes(f, 2), i in axes(f, 1)
        f[i, j, k] = h[i] * (λi / π) * exp(-λi * v[i, j, k]^2) * exp(-λi * w[i, j, k]^2)
    end

    return f

end

"""
$(SIGNATURES)
"""
full_distribution(
    h::X,
    b::X,
    u::Y,
    weights::Y,
    v::Z,
    w::Z,
    prim::A,
    γ = 5 / 3,
) where {X<:AA{<:FN,1},Y<:AA{<:FN,1},Z<:AA{<:FN,3},A<:AA{<:Real,1}} =
    full_distribution(h, b, u, weights, v, w, prim[1], γ)

"""
$(SIGNATURES)

Shift distribution function by external force
"""
function shift_pdf!(f::AV{T}, a, du, dt) where {T<:Real}

    q0 = eachindex(f) |> first # for OffsetArray
    q1 = eachindex(f) |> last

    if a > 0
        shift = Int(floor(a * dt / du)) # only for uniform velocity grid
        for k = q1:-1:q0+shift
            f[k] = f[k-shift]
        end
        for k = q0:shift+q0-1
            f[k] = 0.0
        end

        for k = q0+1:q1
            f[k] += (dt * a - du * shift) * (f[k-1] - f[k]) / du
        end
    else
        shift = Int(floor(-a * dt / du))
        for k = q0:q1-shift
            f[k] = f[k+shift]
        end
        for k = q1-shift+1:q1
            f[k] = 0.0
        end

        for k = q0:q1-1
            f[k] += (dt * a + du * shift) * (f[k] - f[k+1]) / du
        end
    end

    f[q0] = f[q0+1]
    f[q1] = f[q1-1]

    return nothing

end

"""
$(SIGNATURES)

Multi-component gas
"""
function shift_pdf!(f::AM{X}, a::AV{Y}, du::AV{Z}, dt) where {X<:Real,Y<:Real,Z<:Real}

    for j in axes(f, 2)
        _f = @view f[:, j]
        shift_pdf!(_f, a[j], du[j], dt)
    end

    return nothing

end

"""
$(SIGNATURES)

Recover discrete Chapman-Enskog expansion
"""
function chapman_enskog(
    u::AV{T1},
    prim::AV{T2},
    a::AV{T3},
    A::AV{T3},
    τ::Real,
) where {T1<:FN,T2<:Real,T3<:Real}

    M = maxwellian(u, prim)
    f = @. M * (
        1 -
        τ * (a[1] * u + a[2] * u^2 + 0.5 * a[3] * u^3 + A[1] + A[2] * u + 0.5 * A[3] * u^2)
    )

    return f

end

"""
$(SIGNATURES)
"""
function chapman_enskog(
    u::AV{T1},
    prim::AV{<:RN},
    sw::AV{<:RN},
    K::Real,
    τ::Real,
) where {T1}

    Mu, Mxi, _, _1 = gauss_moments(prim, K)
    a = pdf_slope(prim, sw, K)
    swt = -prim[1] .* moments_conserve_slope(a, Mu, Mxi, 1)
    A = pdf_slope(prim, swt, K)

    return chapman_enskog(u, prim, a, A, τ)

end

"""
$(SIGNATURES)
"""
function chapman_enskog(
    u::AA{T1},
    v::AA{T1},
    prim::AV{T2},
    a::AV{T3},
    b::AV{T3},
    A::AV{T3},
    τ::Real,
) where {T1<:FN,T2<:Real,T3<:Real}

    M = maxwellian(u, v, prim)
    f = @. M * (
        1.0 - τ * (a[1] * u + a[2] * u^2 + a[3] * u * v + 0.5 * a[4] * u * (u^2 + v^2)) -
        τ * (b[1] * v + b[2] * u * v + b[3] * v^2 + 0.5 * b[4] * v * (u^2 + v^2)) -
        τ * (A[1] + A[2] * u + A[3] * v + 0.5 * A[4] * (u^2 + v^2))
    )

    return f

end

"""
$(SIGNATURES)
"""
function chapman_enskog(
    u::AA{T1},
    v::AA{T1},
    prim::AV{<:RN},
    swx::AV{<:RN},
    swy::AV{<:RN},
    K::Real,
    τ::Real,
) where {T1}

    Mu, Mv, Mxi, _, _1 = gauss_moments(prim, K)
    a = pdf_slope(prim, swx, K)
    b = pdf_slope(prim, swy, K)
    sw =
        -prim[1] .* (
            moments_conserve_slope(a, Mu, Mv, Mxi, 1, 0) .+
            moments_conserve_slope(b, Mu, Mv, Mxi, 0, 1)
        )
    A = pdf_slope(prim, sw, K)

    return chapman_enskog(u, v, prim, a, b, A, τ)

end

"""
$(SIGNATURES)
"""
function chapman_enskog(
    u::AA{T1},
    v::AA{T1},
    w::AA{T1},
    prim::AV{T2},
    a::AV{T3},
    b::AV{T3},
    c::AV{T3},
    A::AV{T3},
    τ::Real,
) where {T1<:FN,T2<:Real,T3<:Real}

    M = maxwellian(u, v, w, prim)
    f = @. M * (
        1.0 -
        τ * (
            a[1] * u +
            a[2] * u^2 +
            a[3] * u * v +
            a[4] * u * w +
            0.5 * a[5] * u * (u^2 + v^2 + w^2)
        ) -
        τ * (
            b[1] * v +
            b[2] * u * v +
            b[3] * v^2 +
            b[4] * w * v +
            0.5 * b[5] * v * (u^2 + v^2 + w^2)
        ) -
        τ * (
            c[1] * w +
            c[2] * u * w +
            c[3] * v * w +
            c[4] * w^2 +
            0.5 * c[5] * w * (u^2 + v^2 + w^2)
        ) -
        τ * (A[1] + A[2] * u + A[3] * v + A[4] * w + 0.5 * A[5] * (u^2 + v^2 + w^2))
    )

    return f

end

"""
$(SIGNATURES)
"""
function chapman_enskog(
    u::AA{T1},
    v::AA{T1},
    w::AA{T1},
    prim::AV{<:RN},
    swx::AV{<:RN},
    swy::AV{<:RN},
    swz::AV{<:RN},
    K::Real,
    τ::Real,
) where {T1}

    Mu, Mv, Mw, _, _1 = gauss_moments(prim, K)
    a = pdf_slope(prim, swx, K)
    b = pdf_slope(prim, swy, K)
    c = pdf_slope(prim, swz, K)
    sw =
        -prim[1] .* (
            moments_conserve_slope(a, Mu, Mv, Mw, 1, 0, 0) .+
            moments_conserve_slope(b, Mu, Mv, Mw, 0, 1, 0) .+
            moments_conserve_slope(c, Mu, Mv, Mw, 0, 0, 1)
        )
    A = pdf_slope(prim, sw, K)

    return chapman_enskog(u, v, w, prim, a, b, c, A, τ)

end


text:
using Gadfly

#plot(layer(x=[0, 1, 1, 0], y=[0, 1, 0, 1], ygroup=["A", "A", "B", "B"],
           #Geom.subplot_grid(Geom.line), Theme(default_color=colorant"red")),
     #layer(x=[0.5, 1.5, 1.5, 0.5], y=[0, 1, 0, 1], ygroup=["A", "A", "B", "B"],
           #Geom.subplot_grid(Geom.line), Theme(default_color=colorant"blue")))

plot(ygroup=["A", "A", "B", "B"],
     Geom.subplot_grid(
        layer(x=[0, 1, 1, 0], y=[0, 1, 0, 1],
              Geom.line, Theme(default_color=colorant"red")),
        layer(x=[0.5, 1.5, 1.5, 0.5], y=[0, 1, 0, 1],
              Geom.line, Theme(default_color=colorant"blue"))))


text:
"""
Struct for Heston Process

		hstProcess=HestonProcess(σ::num1,σ_zero::num2,λ::num3,κ::num4,ρ::num5,θ::num6) where {num1,num2,num3,num4,num5,num6 <: Number}

Where:\n
		σ	=	volatility of volatility of the process.
		σ_zero	= initial volatility of the process.
		λ	=	??volatility of the process.
		κ	=	??volatility of the process.
		ρ	=	??volatility of the process.
		θ	=	??drift of the process.
"""
mutable struct HestonProcess{num <: Number, num1 <: Number, num2 <: Number, num3 <: Number, num4 <: Number, num5 <: Number, abstrUnderlying <: AbstractUnderlying, numtype <: Number} <: ItoProcess{numtype}
    σ::num
    σ_zero::num1
    λ::num2
    κ::num3
    ρ::num4
    θ::num5
    underlying::abstrUnderlying
    function HestonProcess(σ::num, σ_zero::num1, λ::num2, κ::num3, ρ::num4, θ::num5, underlying::abstrUnderlying) where {num <: Number, num1 <: Number, num2 <: Number, num3 <: Number, num4 <: Number, num5 <: Number, abstrUnderlying <: AbstractUnderlying}
        if σ_zero <= 0.0
            error("initial volatility must be positive")
        elseif σ <= 0.0
            error("volatility of volatility must be positive")
        elseif abs(κ + λ) <= 1e-14
            error("unfeasible parameters")
        elseif !(-1.0 <= ρ <= 1.0)
            error("ρ must be a correlation")
        else
            zero_typed = zero(num) + zero(num1) + zero(num2) + zero(num3) + zero(num4) + zero(num5)
            return new{num, num1, num2, num3, num4, num5, abstrUnderlying, typeof(zero_typed)}(σ, σ_zero, λ, κ, ρ, θ, underlying)
        end
    end
end

export HestonProcess;

function simulate!(X, mcProcess::HestonProcess, rfCurve::ZeroRate, mcBaseData::MonteCarloConfiguration{type1, type2, type3, SerialMode, type5}, T::numb) where {numb <: Number, type1 <: Integer, type2 <: Integer, type3 <: StandardMC, type5 <: Random.AbstractRNG}
    r = rfCurve.r
    S0 = mcProcess.underlying.S0
    d = dividend(mcProcess)
    Nsim = mcBaseData.Nsim
    Nstep = mcBaseData.Nstep
    σ = mcProcess.σ
    σ_zero = mcProcess.σ_zero
    λ1 = mcProcess.λ
    κ = mcProcess.κ
    ρ = mcProcess.ρ
    θ = mcProcess.θ
    @assert T > 0.0

    ####Simulation
    ## Simulate
    κ_s = κ + λ1
    θ_s = κ * θ / κ_s

    dt = T / Nstep
    isDualZero = T * r * σ_zero * θ_s * κ_s * σ * ρ * 0.0 * S0
    #X=Matrix{typeof(isDualZero)}(undef,Nsim,Nstep+1);
    view(X, :, 1) .= isDualZero
    for i = 1:Nsim
        v_m = σ_zero^2
        for j = 1:Nstep
            e1 = randn(mcBaseData.rng)
            e2 = e1 * ρ + randn(mcBaseData.rng) * sqrt(1 - ρ * ρ)
            @views X[i, j+1] = X[i, j] + ((r - d) - 0.5 * v_m) * dt + sqrt(v_m) * sqrt(dt) * e1
            v_m += κ_s * (θ_s - v_m) * dt + σ * sqrt(v_m) * sqrt(dt) * e2
            #when v_m = 0.0, the derivative becomes NaN
            v_m = max(v_m, isDualZero)
        end
    end
    ## Conclude
    X .= S0 .* exp.(X)
    return
end

function simulate!(X, mcProcess::HestonProcess, rfCurve::ZeroRate, mcBaseData::MonteCarloConfiguration{type1, type2, type3, SerialMode, type5}, T::numb) where {numb <: Number, type1 <: Integer, type2 <: Integer, type3 <: AntitheticMC, type5 <: Random.AbstractRNG}
    r = rfCurve.r
    S0 = mcProcess.underlying.S0
    d = dividend(mcProcess)
    Nsim = mcBaseData.Nsim
    Nstep = mcBaseData.Nstep
    σ = mcProcess.σ
    σ_zero = mcProcess.σ_zero
    λ1 = mcProcess.λ
    κ = mcProcess.κ
    ρ = mcProcess.ρ
    θ = mcProcess.θ
    @assert T > 0.0

    ####Simulation
    ## Simulate
    κ_s = κ + λ1
    θ_s = κ * θ / κ_s

    dt = T / Nstep
    isDualZero = T * r * σ_zero * κ * θ * λ1 * σ * ρ * 0.0
    #X=Matrix{typeof(isDualZero)}(undef,Nsim,Nstep+1);
    view(X, :, 1) .= isDualZero
    for i = 1:div(Nsim, 2)
        v_m_1 = σ_zero^2
        v_m_2 = σ_zero^2
        for j = 1:Nstep
            e1 = randn(mcBaseData.rng)
            e2 = -(e1 * ρ + randn(mcBaseData.rng) * sqrt(1 - ρ * ρ))
            @views X[2*i-1, j+1] = X[2*i-1, j] + ((r - d) - 0.5 * v_m_1) * dt + sqrt(v_m_1) * sqrt(dt) * e1
            @views X[2*i, j+1] = X[2*i, j] + ((r - d) - 0.5 * v_m_2) * dt + sqrt(v_m_2) * sqrt(dt) * (-e1)
            v_m_1 += κ_s * (θ_s - v_m_1) * dt + σ * sqrt(v_m_1) * sqrt(dt) * e2
            v_m_2 += κ_s * (θ_s - v_m_2) * dt + σ * sqrt(v_m_2) * sqrt(dt) * (-e2)
            v_m_1 = max(v_m_1, isDualZero)
            v_m_2 = max(v_m_2, isDualZero)
        end
    end
    ## Conclude
    X .= S0 .* exp.(X)
    return
end


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

Reads (length,azimuth) pairs from file and plot a windwiggle diagram.

Full option list at [`pswiggle`](http://gmt.soest.hawaii.edu/doc/latest/pswiggle.html)

Parameters
----------

- $(GMT.opt_R)
- **Z** : **scale** : -- Number or Str --

    Gives anomaly scale in data-units/distance-unit.
	[`-Z`](http://gmt.soest.hawaii.edu/doc/latest/pswiggle.html#z)
- **A** : **azimuth** : -- Str or number --

    Sets the preferred positive azimuth. Positive wiggles will “gravitate” towards that direction.
	[`-A`](http://gmt.soest.hawaii.edu/doc/latest/pswiggle.html#a)
- $(GMT.opt_B)
- **C** : **center** : -- Number --

    Subtract center from the data set before plotting [0].
	[`-C`](http://gmt.soest.hawaii.edu/doc/latest/pswiggle.html#c)
- **D** : **scale_bar** : -- Str --

    Defines the reference point on the map for the vertical scale bar using one of four coordinate systems.
	[`-D`](http://gmt.soest.hawaii.edu/doc/latest/pswiggle.html#d)
- **F** : **bar_rectangle** : -- Str --

    Without further options, draws a rectangular border around the vertical scale bar.
	[`-F`](http://gmt.soest.hawaii.edu/doc/latest/pswiggle.html#f)
- **G** : **fill** : -- Number or Str --

    Set fill shade, color or pattern for positive and/or negative wiggles [Default is no fill].
	[`-G`](http://gmt.soest.hawaii.edu/doc/latest/pswiggle.html#g)
- **I** : **fixed_azim** : -- Number --

    Set a fixed azimuth projection for wiggles [Default uses track azimuth, but see -A].
	[`-I`](http://gmt.soest.hawaii.edu/doc/latest/pswiggle.html#i)
- $(GMT.opt_J)
- $(GMT.opt_Jz)
- $(GMT.opt_P)
- **T** : **pen** : -- Number or Str or Tuple or [] --

    Draw track [Default is no track]. Append pen attributes to use [Defaults: width = 0.25p, color =
    black, style = solid].
	[`-T`](http://gmt.soest.hawaii.edu/doc/latest/pswiggle.html#t)
- **W** : **pen** : -- Number or Str or tuple or [] --

    Specify outline pen attributes [Default is no outline].
	[`-W`](http://gmt.soest.hawaii.edu/doc/latest/pswiggle.html#w)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_swap_xy)
"""
function wiggle(cmd0::String="", arg1=[]; K=false, O=false, first=true, kwargs...)

	length(kwargs) == 0 && return monolitic("pswiggle", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)
	output, opt_T, fname_ext = fname_out(d)		# OUTPUT may have been an extension only

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX12c/12c")
	cmd, = parse_JZ(cmd, d)
	cmd  = parse_UVXY(cmd, d)
	cmd, opt_bi = parse_bi(cmd, d)
	cmd, opt_di = parse_di(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_g(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, opt_i = parse_i(cmd, d)
	cmd, = parse_p(cmd, d)
	cmd, = parse_t(cmd, d)
	cmd, = parse_swap_xy(cmd, d)
	cmd = parse_params(cmd, d)

	cmd, K, O, opt_B = set_KO(cmd, opt_B, first, K, O)		# Set the K O dance

	# If file name sent in, read it and compute a tight -R if this was not provided
	cmd, arg1, opt_R, = read_data(d, cmd0, cmd, arg1, opt_R, opt_i, opt_bi, opt_di)

	cmd = add_opt(cmd, 'A', d, [:A :azimuth])
	cmd = add_opt(cmd, 'C', d, [:C :center])
	cmd = add_opt(cmd, 'D', d, [:D :scale_bar])
	cmd = add_opt(cmd, 'F', d, [:F :bar_rectangle])
	cmd = add_opt(cmd, 'G', d, [:G :fill])
	cmd = add_opt(cmd, 'I', d, [:I :fixed_azim])
	cmd = cmd * opt_pen(d, 'T', [:T :pen])
	cmd = cmd * opt_pen(d, 'W', [:W :pen])
	cmd = add_opt(cmd, 'Z', d, [:Z :scale])

	cmd = finish_PS(d, cmd, output, K, O)
	return finish_PS_module(d, cmd, "", output, fname_ext, opt_T, K, "pswiggle", arg1)
end

# ---------------------------------------------------------------------------------------------------
wiggle!(cmd0::String="", arg1=[]; K=true, O=true,  first=false, kw...) =
	wiggle(cmd0, arg1; K=K, O=O,  first=first, kw...)

wiggle(arg1=[], cmd0::String=""; K=false, O=false,  first=true, kw...) =
	wiggle(cmd0, arg1; K=K, O=O,  first=first, kw...)

wiggle!(arg1=[], cmd0::String=""; K=true, O=true,  first=false, kw...) =
	wiggle(cmd0, arg1; K=K, O=O,  first=first, kw...)

pswiggle  = wiggle			# Alias
pswiggle! = wiggle!			# Alias


text:
module MergingModifiers

println("\nDemonstrating merging modifiers")

using Modia

# Desired:
#   using ModiaMath: plot
#
# In order that these packages need not to be defined in the user environment, they are included via Modia:
using Modia.ModiaMath: plot

@model M begin
  x = Float(start=1.0, info="state")
@equations begin
    der(x) + 2x = 0
  end
end

result = simulate(M, 1, logTranslation=true, removeSingularities=false)
plot(result, ("x"), heading="M", figure=1)

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

@model MInstance begin
  @extends M(x=Var(start=2.0))  # info forgotten
end

result = simulate(MInstance, 1, logTranslation=true, removeSingularities=false)
plot(result, ("x"), heading="MInstance", figure=2)

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

@model MInstance2 begin
  @extends M(y=Var(start=2.0)) # Should give error message for missing y
end

result = simulate(MInstance2, 1, logTranslation=true, removeSingularities=false)
plot(result, ("x"), heading="MInstance2", figure=3)

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

#=
merge(c, e, v) = begin @show c; dump(c); @show e; dump(e.initializers[1].name); dump(e.initializers[1].fdef.args[2].args[2].args[2]); @show v return c end

@model MInstance3 begin
# @extends M(x(start=2.0))  # UndefVarError: x not defined
# @extends M(x=Var(start=2.0))
  @extends M(x=merge(Var(start=2.0), M, :x))
end

result = simulate(MInstance3, 1, logTranslation=true, removeSingularities=false)
plot(result, ("x"), heading="MInstance3", figure=4)
=#

end


text:
# vim:set ft=julia ts=4 sw=4 sts=4 autoindent:

# Run the test suite.
#
# Author:   Pontus Stenetorp    <pontus stenetorp se>
# Version:  2014-12-23

using Hinton

function runtests()
    srand(0x4711)
    w = eye(7, 17) + randn(7, 17)
    open("/dev/null", "w") do f
        println(f, hintontxt(w))
        println(f, hintonvec(w))
    end
end
runtests()


text:
@testset "Unconstrained" begin
  @testset "Basic" begin
    # These basic tests are simple ways of checking that your solver isn't breaking.
    nlp = ADNLPModel(
      x -> (x[1] - 1)^2 + (x[2] - 2)^2 / 4,
      zeros(2)
    )
    output = with_logger(NullLogger()) do
      uncsolver(nlp)
    end
    @test isapprox(output.solution, [1.0; 2.0], rtol=1e-4)
    @test output.objective < 1e-4
    @test output.dual_feas < 1e-4
    @test output.status == :first_order
  end

  @testset "Failures" begin
    # This checks that your solver is failing accordingly
    nlp = ADNLPModel(
      x -> (x[1] - 1)^2 + (x[2] - 2)^2 / 4,
      zeros(2)
    )
    output = with_logger(NullLogger()) do
      uncsolver(nlp, max_eval = 1)
    end
    @test output.status == :max_eval
    output = with_logger(NullLogger()) do
      uncsolver(nlp, max_iter = 1)
    end
    @test output.status == :max_iter
    output = with_logger(NullLogger()) do
      uncsolver(nlp, max_time = 1e-8)
    end
    @test output.status == :max_time
  end
end


text:
module DirectDiffTest

using Test, Random, LinearAlgebra, SparseArrays
using Distributions
import CovSel, ProximalBase, CoordinateDescent

Random.seed!(123)

function genData(p)
  Sigmax = Matrix(1.0I, p, p)
  Sigmay = zeros(p,p)
  rho = 0.7
  for i=1:p
    for j=1:p
      Sigmay[i,j]=rho^abs(i-j)
    end
  end
  sqmy = sqrt(Sigmay)
  n = 1000
  X = randn(n,p)
  Y = randn(n,p) * sqmy;
  hSx = cov(X)
  hSy = cov(Y)
  hSx, hSy
end

@testset "direct_difference_estimation" begin

  @testset "small" begin
    # simple model that differes in only one edge
    p = 10
    hSx, hSy = genData(p)

    tmp = ones(p,p)
    for i=1:5

      # λ = 0.2
      λ = rand(Uniform(0.05, 0.3))
      g = ProximalBase.ProxL1(λ)

      solShoot = CovSel.Alt.differencePrecisionNaive(hSx, hSy, λ, tmp)
      solShoot1 = CovSel.differencePrecisionActiveShooting(Symmetric(hSx), Symmetric(hSy), g)

      @test convert(Matrix, solShoot1) ≈ convert(Matrix, solShoot) atol=1e-5
    end
  end

  @testset "loss" begin
    for rep=1:5
      p = 10
      hSx, hSy = genData(p)

      A = sprand(p, p, 0.15)
      A = (A + A') / 2.
      As = ProximalBase.SymmetricSparseIterate(A)

      @test CovSel.diffLoss(Symmetric(hSx), As, Symmetric(hSy), 2)  ≈ norm(hSx*A*hSy, 2) atol=1e-12
      @test CovSel.diffLoss(Symmetric(hSx), As, Symmetric(hSy), Inf) ≈ norm(hSx*A*hSy, Inf) atol=1e-12
      @test_throws ArgumentError CovSel.diffLoss(Symmetric(hSx), As, Symmetric(hSy), 3)
    end
  end

  @testset "refit" begin
      for rep=1:50
        p = 3
        hSx, hSy = genData(p)

        S = sprand(p, p, 0.2)
        S = S + S' / 2
        indS = findall(x->x!=0, S)
        ilin = LinearIndices(S)[indS]

        opt = CoordinateDescent.CDOptions(;maxIter=5000, optTol=1e-12, randomize=true)
        x = CovSel.differencePrecisionRefit(Symmetric(hSx), Symmetric(hSy), indS)

        x1 = spzeros(p, p)
        x1[ilin] = (0.5 * (kron(hSy, hSx)[ilin, ilin] + kron(hSx, hSy)[ilin, ilin])) \ (hSy - hSx)[ilin]
        @test convert(Matrix, x1) ≈ convert(Matrix, x) atol = 1e-4
      end
  end

end

end


text:
# test differentiable MPC policy on cartpole swingup
include_model("cartpole")
n = model.n
m = model.m
d = model.d

# fast dynamics
@variables x1[1:n], x2[1:n], u1[1:m], w1[1:d], hs[1:1], ts[1:1]

function fd(model::Cartpole{Midpoint, FixedTime}, x⁺, x, u, w, h, t)
    x⁺ - (x + h[1] * f(model, 0.5 * (x + x⁺), u, w))
end

_dynamics = fd(model, x2, x1, u1, w1, hs, ts)
dynamics = eval(build_function(_dynamics, x2, x1, u1, w1, hs, ts)[2])
∇dynamics_x2 = eval(build_function(Symbolics.jacobian(_dynamics, x2), x2, x1, u1, w1, hs, ts)[1])
∇dynamics_x1 = eval(build_function(Symbolics.jacobian(_dynamics, x1), x2, x1, u1, w1, hs, ts)[1])
∇dynamics_u1 = eval(build_function(Symbolics.jacobian(_dynamics, u1), x2, x1, u1, w1, hs, ts)[1])

# fast objective
function objective(x, u, t)
    J = 0.0
    if t < 10
        J += transpose(x) * x
        J += transpose(u) * u
    else
        J += 10.0 * transpose(x) * x
    end
    return J
end

_obj = objective(x1, u1, t)
obj = eval(build_function)
function objective(x)
    J = 0.0
    J += transpose(x) * x
    return J
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using DifferentialEquations
using Plots
function f(du,u,p,t)
    du[1] = -u[1]
end
u0 = [10.0]
const V = 1
prob = ODEProblem(f,u0,(0.0,10.0))

condition(u,t,integrator) = t==4
affect!(integrator) = integrator.u[1] += 10
cb = DiscreteCallback(condition,affect!)

sol = solve(prob,Tsit5(),callback=cb,tstops=[4.0])
plot(sol)

# control problem
function f(du,u,p,t)
    du[1] = -0.5*u[1] + p
    du[2] = -0.5*u[2]
end

v

u0 = [10.0;10.0]
p = 0.0
prob = ODEProblem(f,u0,(0.0,10.0),p)

const tstop = [5.;8.]
sol = solve(prob,Tsit5(),callback = cbs, tstops=tstop)
plot(sol)

#### vector callback

function f(du,u,p,t)
  du[1] = u[2]
  du[2] = -p
  du[3] = u[4]
  du[4] = 0.0
end

function condition(out,u,t,integrator) # Event when event_f(u,t) == 0
  out[1] = u[1]
  out[2] = (u[3] - 10.0)u[3]
end

function affect!(integrator, idx)
  if idx == 1
    integrator.u[2] = -0.9integrator.u[2]
  elseif idx == 2
    integrator.u[4] = -0.9integrator.u[4]
  end
end

cb = VectorContinuousCallback(condition,affect!,2)

u0 = [50.0,0.0,0.0,2.0]
tspan = (0.0,15.0)
p = 9.8
prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob,Tsit5(),callback=cb,dt=1e-3,adaptive=false)
plot(sol, vars = (3,1))

#### my testing
function dynamics!(du,u,p,t)
  α, β, γ, δ = p
  du[1] = (α - β * u[2]) * u[1]
  du[2] = (-γ + δ * u[1]) * u[2]
end

function condition(out, u , t, integrator)
  out[1] = 20-u[1]
  out[2] = 20-u[2]
end

function affect!(integrator, idx)
  if idx == 1
    integrator.u[1] = 0
  end
  if idx == 2
    integrator.u[2] = 0
  end
end

cb = VectorContinuousCallback(condition,affect!,2)

u0 = [100, 100]
p = [0.55, 0.028, 0.84, 0.026]
tspan = (0.0,20.0)
prob = ODEProblem{true}(dynamics!,u0,tspan,p)
sol = solve(prob,Tsit5(), callback = cb)
plot(sol)
@show sol[end]
print(sol)
minimum(sol)

# vector callback with just parameters
function dynamics!(du,u,p,t)
  α, β, γ, δ = p
  du[1] = (α - β * u[2]) * u[1]
  du[2] = (-γ + δ * u[1]) * u[2]
end

const stop1 = [1.]
const stop2 = [3.]

function condition(u,t,integrator)
  t in stop1
end
function condition2(u,t,integrator)
  t in stop2
end
function affect!(integrator)
  integrator.p[4] += 4
end
function affect2!(integrator)
  integrator.p[4] -= 4
end

save_positions = (true,true)
cb = DiscreteCallback(condition, affect!, save_positions=save_positions)
save_positions = (false,true)
cb2 = DiscreteCallback(condition2, affect2!, save_positions=save_positions)

cbs = CallbackSet(cb,cb2)
# set constant stop values
const stop1 = [1.0]
const stop2 = [3.0]

function time_condition1(out,u,t,integrator)
  t in stop1
end
function time_condition2(out,u,t,integrator)
  t in stop2
end

function time_affect1!(integrator, event_idx)
    integrator.p[4] += 4
end
function time_affect2!(integrator)
    integrator.p[4] -= 4
end

save_positions = (true,true)
raise_cb = DiscreteCallback(time_condition1,
                            time_affect1!,
                            save_positions=save_positions)

save_positions = (false,true)
lower_cb = DiscreteCallback(time_condition2,
                            time_affect2!,
                            save_positions=save_positions)

function pop_condition(out,u,t,integrator)
  out[1] = 20-u[1]
  out[2] = 20-u[2]
end

function pop_condition_test(out,u,t,integrator)
  for i in 1:2
    out[i] = 1e-2 - u[i]
  end
end

  out[1] = 20-u[1]
  out[2] = 20-u[2]
end
function pop_affect!(integrator, event_idx)
  if event_idx == 1
    integrator.u[1] = 0
  end
  if event_idx == 2
    integrator.u[2] = 0
  end
end
function pop_affect_test!(integrator, event_idx)
  for i in 1:2
    if event_idx == i
      integrator.u[i] = 0
    end
  end
end

cont_cb1 = VectorContinuousCallback(pop_condition_test,pop_affect_test!,2,save_positions=(true,true))

cb_set = CallbackSet(cb, cb2, cont_cb)

u0 = [100, 100]
p = [0.55, 0.028, 0.84, 0.026]
tspan = (0.0,20.0)
const tstops = [1.;3.]
prob = ODEProblem{true}(dynamics!,u0,tspan,p)
sol = solve(prob,Tsit5(), callback = cont_cb1)
plot(sol, legend = false)

@show sol[end]
print(sol)
minimum(sol)

## calback set with both parameter change and zero condition

function dynamics!(du,u,p,t)
  α, β, γ, δ = p
  du[1] = (α - β * u[2]) * u[1]
  du[2] = (-γ + δ * u[1]) * u[2]
end

function u1_zero_condition!(u,t,integrator) # Event when event_f(u,t) == 0
  0.5-u[1]
end
function u2_zero_condition!(u,t,integrator) # Event when event_f(u,t) == 0
  0.5-u[2]
end

function u1_zero_affect!(integrator)
  integrator.u[1] = 0
end
function u2_zero_affect!(integrator)
  integrator.u[2] = 0
end

# set constant stop values
const stop1 = [5.0]
const stop2 = [8.0]

function inc_mort_cond(u,t,integrator)
  t in stop1
end

# decrease mortality (stop pulse)
function dec_mort_cond(u,t,integrator)
  t in stop2
end

# apply an affect when condition is reached
function affect_inc!(integrator)
  integrator.p[4] += 10
end

function affect_dec!(integrator)
  integrator.p[4] -= 10
end

save_positions = (true, true)
cb_inc = ContinuousCallback(inc_mort_cond, affect_inc!, save_positions=save_positions)
save_positions = (false, true)
cb_dec = ContinuousCallback(dec_mort_cond, affect_inc!, save_positions=save_positions)

#u1_zero_event = ContinuousCallback(u1_zero_condition!,u1_zero_affect!)
#u2_zero_event = ContinuousCallback(u2_zero_condition!,u2_zero_affect!)

cb_set = CallbackSet(cb_inc, cb_dec)

u0 = [100, 100]
p = [0.55, 0.028, 0.84, 0.026]
tspan = (0.0,20.0)
prob = ODEProblem{true}(dynamics!,u0,tspan,p)
const tstop = [5.;8.]
sol = solve(prob,Tsit5(), callback = cb_set, tstops=tstop)
plot(sol)
@show sol[end]
print(sol)
minimum(sol)

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
# script to compute the gray-scale color for the robsutness figure

# Global bounds -
BIG = 2.5
SMALL = 1e-6

# load the robsutness data -
data_array = readdlm("./robustness/alpha_Z1_99.dat")

# transform -
(number_of_samples,number_of_cases) = size(data_array)
gray_scale_array = zeros(number_of_samples,number_of_cases)
for sample_index = 1:number_of_samples

  for case_index = 1:number_of_cases

    value = abs(data_array[sample_index,case_index])
    if (value>BIG)
      value = BIG
    end

    gray_scale_array[sample_index,case_index] = round(255*(1 - (value - SMALL)/(BIG - SMALL)))
  end
end


text:
#!/usr/bin/env julia

using SpectralElements
using LinearAlgebra,Plots, UnPack
#----------------------------------#
function caseSetup!(dfn::Diffusion)

    kx=2.0
    ky=2.0
    kt=2.0
    ν =1.0

    function utrue(x,y,t)
        ut = @. sin(kx*pi*x)*sin.(ky*pi*y)*cos(kt*pi*t)
        return ut
    end

    function setIC!(u,x,y,t)
        u .= utrue(x,y,t)
        return
    end

    function setBC!(ub,x,y,t)
        ub .= @. 0+0*x
        return
    end

    function setForcing!(f,x,y,t)
        ut = utrue(x,y,t)
        f  .= @. ut*((kx^2+ky^2)*pi^2*ν)
        f .+= @. - sin(kx*pi*x)*sin(ky*pi*y)*sin(kt*pi*t)*(kt*pi)
        return
    end

    function setVisc!(ν,x,y,t)
        ν .= @. 1+0*x
        return
    end

    function callback!(dfn::Diffusion)
        @unpack fld,msh = dfn
        @unpack time, istep = dfn.tstep

        ut = utrue(msh.x,msh.y,time[1])
        u  = fld.u
        er = norm(ut-u,Inf)
        println("Step $(istep[1]), Time=$(time[1]), er=$er")
        plt = meshplt(u,msh)
        plt = plot!(zlims=(-1,1))
        display(plt)
        #frame(anim)
        return
    end

    return callback! ,setIC! ,setBC! ,setForcing! ,setVisc!
end

#----------------------------------#
Ex = 5; nr1 = 8;
Ey = 5; ns1 = 8;

m1 = Mesh(nr1,ns1,Ex,Ey)
bc = ['D','D','D','D']
Temperature = Diffusion(bc,m1,Tf=1.0,dt=0.01)

#anim = Animation()
simulate!(Temperature,caseSetup!(Temperature)...)

#gif(anim,"diffusion.gif",fps=30)
#----------------------------------#


text:
#
# Function that writes the gmd3D to a pdb file
#

# If called already with the grid of the gmd3D computed

function gmd3D_write(mysim :: Simulation, grid :: Vector{GMD3DGrid}, output :: String; scale=nothing)

  n = length(grid)

  if scale == nothing
    minrho = grid[1].rho
    maxrho = grid[1].rho
    for i in 2:n
      minrho = min(minrho,grid[i].rho)
      maxrho = max(maxrho,grid[i].rho)
    end
  else
    minrho = scale[1]
    maxrho = scale[2]
  end

  file = open(output,"w")
  println(file,"REMARK  3D representation of MDDF density ")
  println(file,"REMARK  ")
  println(file,"REMARK  Occupancy column contains distance to solute. ")
  println(file,"REMARK  B-factor column contains the scaled density, such that 99.9 is the maximum ")
  println(file,"REMARK  density obsered and 0. is the minimum density observed. The actual limits are: ")
  println(file,"REMARK  Minimum density: $minrho")
  println(file,"REMARK  Maximum density: $maxrho")
  for i in 1:n
    iat = grid[i].atom
    name = mysim.atoms[iat].name
    resname = mysim.atoms[iat].resname
    chain = mysim.atoms[iat].chain
    resnum = mysim.atoms[iat].resid
    x = grid[i].x[1]
    y = grid[i].x[2]
    z = grid[i].x[3]
    b = 99.9*(grid[i].rho - minrho)/(maxrho-minrho)
    occup = grid[i].dmin
    model = 0
    atom = PDBTools.Atom(i,name,resname,chain,resnum,x,y,z,b,occup,model)
    println(file,PDBTools.write_atom(atom))
  end
  close(file)

end

# If called directly with the solute selection and gmd_solute file name

function gmd3D_write(mysim :: Simulation,
                     solute :: Vector{Int64},
                     gmd_solute_file :: String,
                     output :: String;
                     dmax = 7.0, dmin = 1.0, step = 0.5, scale=nothing)

  grid = gmd3D(mysim, solute, gmd_solute_file,
               dmax = 7.0, dmin = 1.0, step = 0.5)

  gmd3D_write(mysim, grid, output, scale=scale)

end


text:
"""
    get_lane_width(veh::Entity, roadway::Roadway)
Returns the width of the lane where `veh` is.
"""
function get_lane_width(veh::Entity, roadway::Roadway)

    lane = get_lane(roadway, veh.state)

    if n_lanes_left(lane, roadway) > 0
        footpoint = get_footpoint(veh)
        lane_left = roadway[LaneTag(lane.tag.segment, lane.tag.lane + 1)]
        return -proj(footpoint, lane_left, roadway).curveproj.t
    else
        return lane.width
    end
end

"""
    get_markerdist_left(veh::Entity, roadway::Roadway)
distance of `veh` to the left marker of the lane
"""
function get_markerdist_left(veh::Entity, roadway::Roadway)
    t = posf(veh.state).t
    lane_width = get_lane_width(veh, roadway)
    return lane_width/2 - t
end

"""
    get_markerdist_left(veh::Entity, roadway::Roadway)
distance of `veh` to the right marker of the lane
"""
function get_markerdist_right(veh::Entity, roadway::Roadway)
    t = posf(veh.state).t
    lane_width = get_lane_width(veh, roadway)
    return lane_width/2 + t
end


text:
export COLORS, MOVE_FORWARD, TURN_LEFT, TURN_RIGHT, UP, DOWN, LEFT, RIGHT, LRUD, EMPTY, WALL, GOAL, GEM, OBSTACLE
export MoveForward, AbstractObject, Empty, Wall, Goal, Door, key, Gem, Obstacle, Agent
export get_color, get_dir, set_dir!

using Crayons
using Colors

const COLORS = (:red, :green, :blue, :magenta, :yellow, :white)

#####
# Objects
#####

abstract type AbstractObject end

Base.show(io::IO, x::AbstractObject) = print(io, Crayon(foreground=get_color(x), reset=true), convert(Char, x))

struct Empty <: AbstractObject end
const EMPTY = Empty()
Base.convert(::Type{Char}, ::Empty) = '⋅'
get_color(::Empty) = :white

struct Wall <: AbstractObject end
const WALL = Wall()
Base.convert(::Type{Char}, ::Wall) = '█'
get_color(::Wall) = :white

struct Goal <: AbstractObject end
const GOAL = Goal()
Base.convert(::Type{Char}, ::Goal) = '♥'
get_color(::Goal) = :red

struct Door{C} <: AbstractObject end
Door(c) = Door{c}()
Base.convert(::Type{Char}, ::Door) = '⩎'
get_color(::Door{C}) where C = C

struct Key{C} <: AbstractObject end
Key(c) = Key{c}()
Base.convert(::Type{Char}, ::Key) = '⚷'
get_color(::Key{C}) where C = C

struct Gem <: AbstractObject end
const GEM = Gem()
Base.convert(::Type{Char}, ::Gem) = '♦'
get_color(::Gem) = :magenta

struct Obstacle <: AbstractObject end
const OBSTACLE = Obstacle()
Base.convert(::Type{Char}, ::Obstacle) = '⊗'
get_color(::Obstacle) = :blue

#####
# Agent
#####

Base.@kwdef mutable struct Agent <: AbstractObject
    color::Symbol=:red
    dir::LRUD
    inventory::Union{Nothing, AbstractObject, Vector}=nothing
end

function Base.convert(::Type{Char}, a::Agent)
    if     a.dir === UP
        '↑'
    elseif a.dir === DOWN
        '↓'
    elseif a.dir === LEFT
        '←'
    elseif a.dir === RIGHT
        '→'
    end
end

get_color(a::Agent) = a.color
get_dir(a::Agent) = a.dir
set_dir!(a::Agent, d) = a.dir = d

struct Transportable end
struct Nontransportable end
const TRANSPORTABLE = Transportable()
const NONTRANSPORTABLE = Nontransportable()
istransportable(::Type{<:Key}) = TRANSPORTABLE
istransportable(::Type{Gem}) = TRANSPORTABLE
istransportable(x::AbstractObject) = istransportable(typeof(x))

(x::Pickup)(a::Agent, o) = x(istransportable(o), a, o)

function (::Pickup)(::Transportable, a::Agent, o::AbstractObject)
    if isnothing(a.inventory)
        a.inventory = o
        true
    elseif a.inventory isa Vector
        i = findfirst(isnothing, a.v)
        if isnothing(i)
            false
        else
            a.inventory[i] = o
            true
        end
    else
        false
    end
end

function (::Drop)(a::Agent)
    if isnothing(a.inventory)
        nothing
    elseif a.inventory isa AbstractObject
        x = a.inventory
        a.inventory = nothing
        x
    elseif a.inventory isa Vector
        i = findlast(x -> x isa AbstractObject, a.inventory)
        if isnothing(i)
            nothing
        else
            x = a.inventory[i]
            a.inventory[i] = nothing
            x
        end
    else
        @error "unknown inventory type $(a.inventory)"
    end
end


text:
#=
diffs.jl
finite differences
2019-03-06 Jeff Fessler, University of Michigan
2020-06 N-D version by Steven Whitaker
=#

export diff_map, diff_forw, diff_adj

using LinearMapsAA: LinearMapAA
using Test: @test, @test_throws

# size after apply `diff` along `dim`
diff_size(N::Dims, dim::Int) = (N[1:dim-1]..., N[dim] - 1, N[dim+1:end]...)

# corresponding length
diff_length(N::Dims, dim::Int) = prod(diff_size(N, dim))

# check argument validity
# dims should be Int or AbstractVector{Int} or Dims
# but using duck typing for simplicity
function diff_check(N::Dims, dims)
    (length(unique(dims)) != length(dims)) &&
        throw(ArgumentError("non-unique dims $dims"))
	map(dim -> diff_check(N, dim), dims) # check each dim
end

function diff_check(N::Dims, dim::Int)
    !(1 <= dim <= length(N)) && throw(ArgumentError("dim range $dim"))
    (N[dim] == 1) && throw(ArgumentError("invalid size $N for dim $dim"))
end

"""
    d = diff_forw(X ; dims = 1:ndims(X))

Finite differences along one or more dimensions of an array,
e.g., for anisotropic TV regularization.

By default performs the same operations as
``d = [(I_{N_d} \\otimes \\cdots \\otimes D_{N_1}); \\dots; (D_{N_d} \\otimes \\cdots \\otimes I_{N_1})] vec(X)``
where ``D_N`` denotes the `N-1 × N` 1D finite difference matrix
and `⊗` denotes the Kronecker product, but does it efficiently
without using `spdiagm` (or any `SparseArrays` function).

Input dimension `N` must exceed `1` for each dimension specified by `dims`.

in
- `X` `N_1 × ... × N_d` array (typically an N-D image).

option
- `dims` dimension(s) for performing finite differences; default `1:ndims(X)`
must have unique elements and be a nonempty subset of `1:ndims(X)`

out
- `d` vector of default length `N_d*...*(N_1-1) + ... + (N_d-1)*...*N_1`
"""
function diff_forw(x::AbstractArray{<:Number,D} ; dims = 1:D) where {D}
    diff_check(size(x), dims)
    return reduce(vcat, vec(diff(x, dims = d)) for d in dims)
end

"""
    Z = diff_adj(dx, N::Dims{D} ; dims = 1:D)

Adjoint of finite differences of arrays along one or more dimensions.
By default performs the same operations as
``vec(Z) = [(I_{N_d} \\otimes \\cdots \\otimes D_{N_1}); \\dots; (D_{N_d} \\otimes \\cdots \\otimes I_{N_1})]' * d``
where `D_N` denotes the `N-1 × N` 1D finite difference matrix
and `⊗` denotes the Kronecker product,
but does it efficiently without using `spdiagm` (or any `SparseArrays` function).

in
- `dx` vector of typical length `N_d*...*(N_1-1) + ... + (N_d-1)*...*N_1`
- `N::Dims` desired output size

option
- `dims` dimension(s) for performing adjoint finite differences; default `1:ndims(X)`

out
- `Z` `N_1 × ... × N_d` array by default

"""
function diff_adj(d::AbstractVector{<:Number}, N::Dims{D} ; dims = 1:D) where {D}

    length(d) != sum(diff_length(N,dim) for dim in dims) && throw("length(d)")

    z = zeros(eltype(d), N)
    for (i, dim) in enumerate(dims)
        if i == 1
            di = @view(d[1:diff_length(N,dim)])
        else
            start = 1 + sum(diff_length(N,n) for n in dims[1:i-1])
            len = diff_length(N,dim)
            di = @view(d[start:start+len-1])
        end
        di = reshape(di, diff_size(N,dim))
        slice1 = selectdim(z, dim, 1)
        slice1 .-= selectdim(di, dim, 1)
        slicen = selectdim(z, dim, 2:N[dim]-1)
        slicen .+= selectdim(di, dim, 1:N[dim]-2) - selectdim(di, dim, 2:N[dim]-1)
        sliceN = selectdim(z, dim, N[dim])
        sliceN .+= selectdim(di, dim, N[dim] - 1)
    end

    return z
end

# backward compatibility, undocumented because deprecated
function diff2d_forw(x::AbstractMatrix{<:Number})
    isinteractive() && @warn("diff2d_forw is deprecated; use diff_forw")
    return diff_forw(x)
end

function diff2d_adj(d::AbstractVector{<:Number}, M::Int, N::Int ; out2d=false)
    isinteractive() && @warn("diff2d_adj is deprecated; use diff_adj")
    tmp = diff_adj(d, (M, N))
    return out2d ? tmp : vec(tmp)
end

function diff2d_map(M::Int, N::Int)
    isinteractive() && @warn("diff2d_map is deprecated; use diff_map")
    return diff_map((M, N))
end

"""
    T = diff_map(N::Dims{D} ; dims = 1:D)

in
- `N::Dims` image size

out
- `T` `LinearMapAA` object for computing finite differences via `T*x`
"""
function diff_map(N::Dims{D} ; dims = 1:D) where {D}
    diff_check(N, dims)
    return LinearMapAA(
        x -> diff_forw(reshape(x, N), dims=dims),
        d -> vec(diff_adj(d, N, dims=dims)),
        (sum(diff_length(N,dim) for dim in dims), prod(N)),
        (name="diff_map", dims=dims),
    )
end

"""
    diff_map(:test)
self test
"""
function diff_map(test::Symbol)
    test != :test && throw(ArgumentError("test $test"))
    @test_throws ArgumentError diff_forw(ones(3), dims=2)
    @test_throws ArgumentError diff_forw(ones(3), dims=[1,2])
    @test_throws ArgumentError diff_forw(ones(1,2,1), dims=(1,2))
    @test diff_forw(ones(2,4,6) ; dims=2) == zeros(2*3*6)

    for N in [(3,), (3,4), (2,3,4), (4,3,2,2)]
        T = diff_map(N)
        @test Matrix(T)' == Matrix(T')
        @test T.name == "diff_map"
        T = diff_map(N, dims=[1])
        @test Matrix(T)' == Matrix(T')
        if length(N) >= 2
            for dims in [2, (2,), [1,2], (1,2)]
                T = diff_map(N, dims=dims)
                @test Matrix(T)' == Matrix(T')
            end
        end
        if length(N) >= 3
            for dims in [3, (1,3), [2,3], (1,2,3)]
                T = diff_map(N, dims=dims)
                @test Matrix(T)' == Matrix(T')
            end
        end
    end

    if true # test old 2D versions
        N = (3,5)
        x = rand(N...)
        @test diff2d_forw(x) == diff_forw(x)
        d = rand(sum(diff_length(N,dim) for dim=1:2))
        @test diff2d_adj(d, N... ; out2d=false) == vec(diff_adj(d, N))
        @test diff2d_adj(d, N... ; out2d=true) == diff_adj(d, N)
        @test Matrix(diff2d_map(N...)) == Matrix(diff_map(N))
    end

    N = (1,2)
    @test_throws ArgumentError T = diff_map(N)
    @test_throws ArgumentError T = diff_map(N, dims=1)
    T = diff_map(N, dims=2)
    @test Matrix(T)' == Matrix(T')
    true
end