text:
module Logger

import Dates

import Dance.Configuration

"""
    log(error)

Logging of timestamp & error message

- If Prod: Write to log file
- Else: Output te REPL
"""
function log(error::String) :: Nothing
    error::String = string(Dates.now(Dates.UTC), ": ", error)

    if Configuration.is_prod()
        log_filepath::String = Configuration.Settings[:log_filename]*".log"

        # Create logfile dir if not found
        if !isfile(log_filepath)
            paths::Array{String} = split(log_filepath, '/')
            mkpath(replace(log_filepath, paths[end] => ""))
        end

        open(log_filepath, "a+") do file
            write(file, error * "\n")
        end
    else
        @error error
    end
end

end


text:
export RedisPipe

struct RedisPipe{T} <: AbstractRedisCollection
    conn::AbstractRedisConnection
    key::String

    RedisPipe{T}(conn, key) where T = if serializeable(T)
        new(conn, key)
    else
        throw(ArgumentError("RedisPipe currently not supports arbitrary element type"))
    end

    RedisPipe{T}(key) where T = RedisPipe{T}(default_connection, key)
end

reply(rp::RedisPipe, x) = throw(ProtocolException("unexpected return value $x"))
reply(rp::RedisPipe, x::Vector) = reply(rp, x[2])
reply(rp::RedisPipe, x::String) = reply(rp, x.data)
reply(rp::RedisPipe, ::Nothing) = throw(TimeoutException("timeout"))
reply(rp::RedisPipe{T}, x::Bytes) where T = deserialize(T, x)

function read(rp::RedisPipe, timeout::Integer=0)
    timeout < 0 && throw(ArgumentError("timeout must be non-negative"))
    res = exec(rp.conn, "blpop", rp.key, timeout)
    reply(rp, res)
end

take!(rp::RedisPipe, timeout::Integer=0) = read(rp, timeout)

function peek(rp::RedisPipe)
    res = exec(rp.conn, "lindex", rp.key, 0)
    reply(rp, @some(res))
end

function write(rp::RedisPipe{T}, value) where T
    exec(rp.conn, "rpush", rp.key, serialize(T, value))
    rp
end

put!(rp::RedisPipe, value) = write(rp, value)

function length(rp::RedisPipe)
    exec(rp.conn, "llen", rp.key)::Int64
end

function isempty(rp::RedisPipe)
    length(rp) == 0
end

ext: jl


text:
#!/usr/bin/env julia
#=
Utility functions
del2z <delta.z@aliyun.com>
=#
module Utils

using ..Ant: Polar, Model
export nothing

end # module


text:
StatsBase.vcov(rm::RegressionModel) = scale(rm,true) * unscaledvcov(rm)

type NonlinearLS{T<:FP} <: RegressionModel # nonlinear least squares fits
    m::NLregModF{T}
    pars::Vector{T}
    incr::Vector{T}
    ch::Base.Cholesky{T}
    rss::T      # residual sum of squares at last successful iteration
    tolsqr::T   # squared tolerance for orthogonality criterion
    minfac::T
    mxiter::Int
    fit::Bool
end
function NonlinearLS{T<:FP}(m::NLregModF{T},init::Vector{T})
    p,n = size(m)
    if isa(m,PLregMod)
        nnl,nl,n = size(m)
        p = nl + nnl
    end
    length(init) == p || error("Dimension mismatch")
    rss = updtmu!(m, init); tg = m.tgrad
    NonlinearLS(m, init, zeros(T,p), cholfact(eye(p)), rss, 1e-8, 0.5^10, 1000, false)
end
NonlinearLS{T<:FP}(m::NLregModF{T}) = NonlinearLS(m, initpars(m))

## returns a copy of the current parameter values
StatsBase.coef(nl::NonlinearLS) = copy(nl.pars)

## returns the coefficient table
function StatsBase.coeftable(nl::NonlinearLS)
    pp = coef(nl); se = stderr(nl); tt = pp ./ se
    CoefTable(hcat(pp, se, tt, ccdf(FDist(1, df_residual(nl)), abs2(tt))),
              ["Estimate","Std.Error","t value", "Pr(>|t|)"],
              map(string,pnames(nl)[1:length(pp)]), 4)
end

StatsBase.deviance(nl::NonlinearLS) = nl.rss
df_residual(nl::NonlinearLS) = ((p,n) = size(nl);n-p)
function gnfit(nl::NonlinearLS,verbose::Bool=false) # Gauss-Newton nonlinear least squares
    if !nl.fit
        m = nl.m; pars = nl.pars; incr = nl.incr; minf = nl.minfac; cvg = 2(tol = nl.tolsqr)
        r = m.resid; tg = m.tgrad; ch = nl.ch; nl.rss = rss = updtmu!(m,pars); UL = ch.UL
        for i in 1:nl.mxiter
            ## Create the Cholesky factor of tg * tg' in place
            _,info = LAPACK.potrf!('U',BLAS.syrk!('U','N',1.,tg,0.,UL))
            info == 0 || error("Singular gradient matrix at pars = $(pars')")
            ## solve in place for the Gauss-Newton increment - done in two stages
            ## to be able to evaluate the orthogonality convergence criterion
            cvg = sumsq(BLAS.trsv!('U','T','N',UL,BLAS.gemv!('N',1.,tg,r,0.,incr)))/rss
            if verbose
                print("Iteration:",lpad(string(i),3),", rss = "); showcompact(rss)
                print(", cvg = "); showcompact(cvg); print(" at [");
                showcompact(pars[1])
                for i in 2:length(pars); print(","); showcompact(pars[i]);end
                println("]")
            end
            BLAS.trsv!('U','N','N',UL,incr)
            if verbose
                print("   Incr: [")
                showcompact(incr[1])
                for i in 2:length(incr); print(","); showcompact(incr[i]);end
                println("]")
            end
            f = 1.
            while true
                f >= minf || error("Failure to reduce rss at $(pars') with incr = $(incr') and minfac = $minf")
                rss = updtmu!(nl.m, pars + f * incr)
                if verbose
                    print("  f = ",f,", rss = "); showcompact(rss); println()
                end
                rss < nl.rss && break
                f *= 0.5                    # step-halving
            end
            cvg < tol && break
            pars += f * incr
            nl.rss = rss
        end
        verbose && println()
        copy!(nl.pars,pars)
        cvg < tol || error("failure to converge in $(nl.mxiter) iterations")
        nl.fit = true
    end
    nl
end
model(nl::NonlinearLS) = nl.m
StatsBase.nobs(nl::NonlinearLS) = size(nl,2)
pnames(nl::NonlinearLS) = pnames(model(nl))

function Base.show{T<:FP}(io::IO, nl::NonlinearLS{T})
    gnfit(nl)
    p,n = size(nl)
    s2 = deviance(nl)/convert(T,n-p)
    println(io, "Nonlinear least squares fit to ", n, " observations")
    ## Add a model or modelformula specification in here
    println(io); show(io, coeftable(nl)); println(io)
    print(io,"Residual sum of squares at estimates: "); showcompact(io,nl.rss); println(io)
    print(io,"Residual standard error = ");showcompact(io,sqrt(s2));
    print(io, " on $(n-p) degrees of freedom")
end
Base.size(nl::NonlinearLS) = size(tgrad(model(nl)))
Base.size(nl::NonlinearLS,args...) = size(tgrad(model(nl)),args...)

type PLinearLS{T<:FP} <: RegressionModel
    m::PLregModF{T}
    qr::QRCompactWY{T}
    pars::Vector{T}
    incr::Vector{T}
    B::Matrix{T}
    rss::T      # residual sum of squares at last successful iteration
    tolsqr::T   # squared tolerance for orthogonality criterion
    minfac::T
    mxiter::Int
    fit::Bool
end
function PLinearLS{T<:FP}(m::PLregModF{T},nlpars::Vector{T})
    nnl,nl,n = size(m)
    length(nlpars) == nnl || throw(DimensionMismatch(""))
    qr = qrfact!(updtMM!(m,nlpars)')
    pars = [vec(qr\model_response(m)),nlpars]
    PLinearLS(m, qr, pars, zeros(T,nnl), Array(T,n,nnl), updtmu!(m,pars),
              convert(T,1e-8), convert(T,0.5^10), 1000, false)
end
PLinearLS(m::PLregModF) = PLinearLS(m,initpars(m))

pnames(pl::PLinearLS) = pnames(pl.m)

function StatsBase.deviance{T<:FP}(pl::PLinearLS{T},nlp::Vector{T})
    m = pl.m; nnl,nl,n = size(m); pars = pl.pars
    copy!(view(pars,nl + (1:nnl)), nlp) # record nl pars
    pl.qr = qr = qrfact!(updtMM!(m,nlp)') # update and decompose lin pars model matrix
    copy!(view(pars,1:nl),qr\model_response(m)) # conditionally optimal linear pars
    pl.rss = updtmu!(m,pars)         # update mu and evaluate deviance
end
StatsBase.deviance(pl::PLinearLS) = pl.rss

function gpinc{T<:FP}(pl::PLinearLS{T})
    m = pl.m; nnl,nl,n = size(m); Aphi = mmjac(m); B = pl.B
    r = residuals(m); mqr = pl.qr
    lin = 1:nl; lpars = pl.pars[lin]
    for k in 1:nnl
        B[:,k] = reshape(Aphi[k,:,:],(nl,n))' * lpars
    end
    LAPACK.gemqrt!('L','T',mqr.factors,mqr.T.data,B)
    for j in 1:nnl, i in 1:nl
        B[i,j] = dot(vec(Aphi[j,i,:]),r)
    end
    BLAS.trsm!('L','U','N','N',1.0,mqr[:R],sub(B,lin,:))
    rhs = mqr[:Q]'*model_response(m); for i in 1:nl rhs[i] = zero(T) end
    st = qrfact(B); sc = (st[:Q]'*rhs)[1:nnl]
    BLAS.trsv!('U','N','N',sub(st.factors,1:nnl,:),copy!(pl.incr,sc))
    sumabs2(sc)/pl.rss
end

function gpfit(pl::PLinearLS,verbose::Bool=false) # Golub-Pereyra variable projection algorithm
    if !pl.fit
        m = pl.m; pars = pl.pars; nnl,nl,n = size(m); nlin = nl + (1:nnl);
        minf = pl.minfac; cvg = 2(tol = pl.tolsqr); oldrss = pl.rss
        for i in 1:pl.mxiter
            cvg = gpinc(pl) # evaluate the increment and convergence criterion
            if verbose
                print("Iteration:",lpad(string(i),3),", rss = "); showcompact(oldrss)
                print(", cvg = "); showcompact(cvg); print(" at [")
                showcompact(pl.pars[1])
                for i in 2:length(pl.pars); print(","); showcompact(pl.pars[i]); end
                println("]")
                print("   Incr: [")
                showcompact(pl.incr[1])
                for i in 2:length(pl.incr); print(","); showcompact(pl.incr[i]); end
                println("]")
            end
            f = 1.
            while true
                f >= minf || error("Failure to reduce rss at $(pars') with incr = $(pl.incr') and minfac = $minf")
                deviance(pl, pl.pars[nlin] + f * pl.incr)
                if verbose
                    print("  f = ",f,", rss = "); showcompact(pl.rss); println()
                end
                pl.rss < oldrss && break
                f *= 0.5                    # step-halving
            end
            cvg < tol && break
            oldrss = pl.rss
        end
        verbose && println()
        cvg < tol || error("failure to converge in $(pl.mxiter) iterations")
        pl.fit = true
    end
    pl
end

## returns a copy of the current parameter values
StatsBase.coef(pl::PLinearLS) = copy(pl.pars)

## returns the coefficient table
function StatsBase.coeftable(pl::PLinearLS)
    pp = coef(pl); se = stderr(pl); tt = pp ./ se
    CoefTable(hcat(pp, se, tt, ccdf(FDist(1, df_residual(pl)), abs2(tt))),
              ["Estimate","Std.Error","t value", "Pr(>|t|)"],
              map(string,pnames(pl)[1:length(pp)]), 4)
end

df_residual(pl::PLinearLS) = nobs(pl) - npars(pl)

StatsBase.fit(m::PLregModF,verbose=false) = gpfit(PLinearLS(m),verbose)

StatsBase.model_response(pl::PLinearLS) = model_response(pl.m)

npars(pl::PLinearLS) = npars(pl.m)

function Base.scale(pl::PLinearLS,sqr=false)
    scsq = deviance(pl)/df_residual(pl)
    sqr ? scsq : sqrt(scsq)
end

StatsBase.residuals(pl::PLinearLS) = residuals(pl.m)

function Base.show{T<:FP}(io::IO, pl::PLinearLS{T})
    gpfit(pl)
    println(io, "Nonlinear least squares fit to $(nobs(pl)) observations")
## Add a model or modelformula specification in here
    println(io); show(io, coeftable(pl)); println(io)
    print(io,"Residual sum of squares at estimates: "); showcompact(io,pl.rss); println(io)
    print(io,"Residual standard error = ");showcompact(io,scale(pl));
    print(io, " on $(df_residual(pl)) degrees of freedom")
end

Base.size(pl::PLinearLS) = size(pl.m.MMD)

Base.size(pl::PLinearLS,args...) = size(pl.m.MMD,args...)

StatsBase.vcov(pl::PLinearLS) = scale(pl,true) * unscaledvcov(pl.m)


text:
using DelimitedFiles, LightGraphs

####################################################################################
#######################  Increment Counter and Get Next Dir  ########################
####################################################################################

function nextRun(dir)
    dir = joinpath(@__DIR__, "./$dir")
    if isfile("$(dir)/counter.txt") == false
        mkpath("$(dir)")
        writedlm("$(dir)/counter.txt", "1")
    end
    c = Int(readdlm("$(dir)/counter.txt")[1])
    writedlm("$(dir)/counter.txt", c + 1)
    newdir = "$(dir)/$(c)"
    mkdir(newdir)
    return newdir
end

####################################################################################
#############################  Read Topology  ######################################
####################################################################################

function readTopology(topology; zeroindex=false, downscale=1, what_to_read="topology.txt")
    dir = joinpath(@__DIR__, "./data/$topology")
    input_topology = readdlm("$dir/$what_to_read", header=true)[1]
    input_nodes = readdlm("$dir/nodes.txt", header=true)[1]

    ignore = ()
    fromNodes = input_topology[:,1]
    toNodes = input_topology[:,2]
    capacity = input_topology[:,3] ./ downscale ./ 1000
    probabilities = input_topology[:,4]
    links = []
    for i in 1:size(fromNodes, 1)
        if (!(fromNodes[i] in ignore) && !(toNodes[i] in ignore))
            push!(links, (Int(fromNodes[i] + zeroindex), Int(toNodes[i] + zeroindex)))
        end
    end
    return links, capacity, probabilities, input_nodes
end

function writeTopology(topology, links, capacity, probabilities)
	dir = joinpath(@__DIR__, "./data/$topology")
	ntopo = open("$dir/topo_n.txt", "w")
	# this header is incorrect; wtf
	print(ntopo, "to_node from_node capacity prob_failure\n")
	for i in 1:length(links)
		print(ntopo, links[i][1], "\t", links[i][2], "\t", capacity[i]*1000, "\t", probabilities[i], "\n")
	end
	flush(ntopo)
	close(ntopo)
	return
end

####################################################################################
#############################  Read Demand  ########################################
####################################################################################

function readDemand(filename, num_nodes, num_demand; scale=1.0, matrix=true, downscale=1, sigfigs=1, zeroindex=false)
    filename = joinpath(@__DIR__, "./data/$filename")
    input_demand = matrix ? ParseMatrix("$(filename).txt", num_nodes, num_demand) : IgnoreCycles(readdlm("$(filename)/$(num_demand).txt", header=true)[1], zeroindex=zeroindex)
    fromNodes = input_demand[:,1]
    toNodes = input_demand[:,2]
    flows = map(tup -> (Int(tup[2]), Int(toNodes[tup[1]])), enumerate(fromNodes))
    # demand = map(i -> round(i / downscale * sigfigs) / sigfigs * scale + 1/sigfigs, input_demand[:,3])
    demand = input_demand[:,3] ./ downscale .* scale ./ 1000
    return demand, flows
end

function IgnoreCycles(demand; zeroindex=false)
    z = [0 0 0;]
    for row in 1:size(demand,1)
        if demand[row, 1] != demand[row, 2]
            z = vcat(z, [demand[row,1] + zeroindex demand[row, 2] + zeroindex demand[row,3];])
            # z = vcat(z, transpose(demand[row,:]))
        end
    end
    return z[2:end,:]
end

####################################################################################
#############################  Get Allocations  ####################################
####################################################################################
function parseYatesSplittingRatios(filename, k, flows; zeroindex=false)
    filename = joinpath(@__DIR__, "./data/$filename")
    f = readdlm(filename)
    a = zeros(size(flows,1), k)
    num_flow = 0
    t = 0
    for row in 1:size(f,1)
        if "->" in f[row,:]
            fromNode = parse(Int, replace(f[row, 1], "h" => "")) + zeroindex
            toNode = parse(Int, replace(f[row, 3], "h" => "")) + zeroindex
            num_flow = findfirst(x -> x == (fromNode, toNode), flows)
            t = 1
        else
            for col in 1:size(f,2)
                if occursin("@", f[row,col])
                    a[num_flow, t] = ceil(f[row,col + 1] * 1000)/1000
                    t += 1
                    break
                end
            end
        end
    end
    return a
end

function parseYatesAllocations(filename, k, demand, flows; zeroindex=false)
    f = joinpath(@__DIR__, "./data/$filename")
    a = zeros(size(flows,1), k)
    num_flow = 0
    t = 0
    for row in 1:size(f,1)
        if "->" in f[row,:]
            fromNode = parse(Int, replace(f[row, 1], "h" => "")) + zeroindex
            toNode = parse(Int, replace(f[row, 3], "h" => "")) + zeroindex
            num_flow = findfirst(x -> x == (fromNode, toNode), flows)
            t = 1
        else
            for col in 1:size(f,2)
                if occursin("@", f[row,col])
                    a[num_flow, t] = ceil(f[row,col + 1] * demand[num_flow] * 1000)/1000
                    t += 1
                    break
                end
            end
        end
    end
    return a
end

####################################################################################
################################  Parse Paths  #####################################
####################################################################################

function parsePaths(filename, links, flows; zeroindex=false)
    filename = joinpath(@__DIR__, "./data/$filename")
    nflows = length(flows)
    x = readdlm(filename)
    T = []
    Tf = [[] for i=1:nflows]
    tf = []
    fromNode = 0
    toNode = 0
    num_flow = 0
    tindex = 1
    max_paths = 0
    paths = Matrix(undef, nflows, nflows)
    for row in 1:size(x,1)
        if "->" in x[row,:]
            max_paths = max(max_paths, size(tf, 1))
            if fromNode != 0 && num_flow != 0
                #add tunnel if not first run
                paths[fromNode, toNode] = tf
                Tf[num_flow] = tf
            end
            if x[row,1] isa Number
                fromNode = x[row, 1] + zeroindex
                toNode = x[row, 3] + zeroindex
            else
                fromNode = parse(Int, replace(x[row, 1], "h" => "")) + zeroindex
                toNode = parse(Int, replace(x[row, 3], "h" => "")) + zeroindex
            end
            num_flow = findfirst(x -> x == (fromNode, toNode), flows)
            if num_flow == nothing
                num_flow = 0
            end
            tf = []
        else
            t = []
            #parse each tunnel into edges
            for col in 2:size(x,2)
                if occursin("]", x[row,col])
                    break
                end
                r = replace(x[row,col][2:end-2], "s" => "")
                stringtup = split(r, ",")
                #find edge in edge matrix
                e = (parse(Int,stringtup[1]) + zeroindex, parse(Int,stringtup[2]) + zeroindex)
                index = findfirst(x -> x == e, links)
                push!(t, index)
            end
            #create new tunnel and add index to that flows tunnels
            if num_flow != 0
                push!(T, t)
                push!(tf, tindex)
                tindex += 1
            end
        end
    end

    #add last tunnel
    if num_flow != 0
        Tf[num_flow] = tf
        paths[fromNode, toNode] = tf
    end
    push!(T, [])

    for f in 1:size(Tf,1)
        for t in 1:max_paths
            try Tf[f][t]
            catch
                push!(Tf[f], tindex)
            end
        end
    end
    return T, Tf, max_paths
end

####################################################################################
################################  Parse Matrix  ####################################
####################################################################################

function ParseMatrix(filename, num_nodes, num_demand)
    ignore = ()
    start_range = 0
    end_range = Inf
    x = readdlm(filename)[num_demand,:]
    m = zeros(length(x)-Int(sqrt(length(x))), 3)
    fromNode = 0
    count = 1
    for i in 0:(num_nodes^2-1)
        toNode = i%num_nodes + 1
        if toNode == 1
            fromNode += 1
        end
        if (fromNode != toNode && i >= start_range && i < end_range && !(fromNode in ignore) && !(toNode in ignore))
            m[count,1] = fromNode
            m[count,2] = toNode
            m[count,3] = x[i+1]
            count += 1
        end
    end
    ret = m
    for row in size(m,1):-1:1
        if (m[row,3] == 0)
            ret = ret[setdiff(1:end, row), :]
        end
    end
    return ret
end

####################################################################################
################################  KSP Tunnels  #####################################
####################################################################################

function getTunnels(nodes, edges, capacity, flows, k; edge_disjoint=false)
    num_edges = length(edges)
    num_nodes = length(nodes)
    num_flows = length(flows)
    graph = LightGraphs.DiGraph(num_nodes)
    distances = Inf*ones(num_nodes, num_nodes)

    for i in 1:num_edges
        LightGraphs.add_edge!(graph, edges[i][1], edges[i][2])
        distances[edges[i][1], edges[i][2]] = 1
        distances[edges[i][2], edges[i][1]] = 1
    end
    T = [[]]
    Tf = []
    ti = 2
    max_k = 1
    for f in 1:num_flows
        tf = []
        curr_k = 0
        state = LightGraphs.yen_k_shortest_paths(graph, flows[f][1], flows[f][2], distances, k)
        paths = state.paths
        edges_used = []
        for i in 1:k
            t = []
            path = i <= size(paths,1) ? paths[i] : []
            for n in 2:size(path,1)
                e = findfirst(x -> x == (path[n-1],path[n]), edges)
                push!(t, e)
            end
            if length(t) == 0 break end
            add = true
            for e in t
                if e in edges_used && edge_disjoint
                    add = false
                    break
                end
            end
            if add
                if edge_disjoint edges_used = vcat(edges_used, t) end
                push!(T, t)
                push!(tf, ti)
                ti += 1
                curr_k += 1
            end
        end
        max_k = max(max_k, curr_k)
        push!(Tf, tf)
    end

    for f in 1:num_flows
        for i in 1:max_k-length(Tf[f])
            push!(Tf[f], 1)
        end
    end
    return T, Tf, max_k, graph
end


text:
"""
    Attention(in => out)
    Attention(qkvlayer)

Self attention layer used by transformer models. Specify the `in` and `out` dimensions,
or directly provide a `qkvlayer` that maps an input the queries, keys, and values.
"""
struct Attention{T}
  qkv::T
end

Attention(dims::Pair{Int, Int}) = Attention(Dense(dims.first, dims.second * 3; bias = false))

@functor Attention

function (attn::Attention)(x::AbstractArray{T, 3}) where T
  q, k, v = chunk(attn.qkv(x), 3, dims = 1)
  scale = convert(T, sqrt(size(q, 1)))
  score = softmax(batched_mul(batched_transpose(q), k) / scale)
  attention = batched_mul(v, score)

  return attention
end

struct MHAttention{S, T}
  heads::S
  projection::T
end

"""
    MHAttention(in, hidden, nheads; dropout = 0.0)

Multi-head self-attention layer used in many vision transformer-like models.

# Arguments
- `in`: Number of dimensions in the input.
- `hidden`: Number of dimensions in the intermediate layer.
- `nheads`: Number of attention heads.
- `dropout`: Dropout rate for the projection layer.
"""
function MHAttention(in, hidden, nheads; dropout = 0.)
  if (nheads == 1 && hidden == in)
    return Attention(in => in)
  end
  inheads, innerheads = chunk(1:in, nheads), chunk(1:hidden, nheads)
  heads = Parallel(vcat, [Attention(length(i) => length(o)) for (i, o) in zip(inheads, innerheads)]...)
  projection = Chain(Dense(hidden, in), Dropout(dropout))

  MHAttention(heads, projection)
end

@functor MHAttention

function (mha::MHAttention)(x)
  nheads = length(mha.heads.layers)
  xhead = chunk(x, nheads, dims = 1)
  return mha.projection(mha.heads(xhead...))
end


text:
mutable struct TauLeapSimulator{M,F} <: AbstractSimulator
  algorithm::M
  next_leap_jumps::Vector{Int}
  next_leap_time::Float64
  t::Float64
  tfinal::Float64
  execute_leap!::F
end

###### constructors #####
function TauLeapSimulator(algorithm::M, number_jumps, execute_leap!::F) where {M,F}
  next_leap_jumps = zeros(Int, number_jumps)

  return TauLeapSimulator{M,F}(algorithm, next_leap_jumps, 0.0, 0.0, 0.0, execute_leap!)
end

##### main simulation routines #####

@inline cumulative_intensity(simulator::TauLeapSimulator) = first(simulator.algorithm.total_rate)

# start a new simulation and initialize the stepper with the next event
@inline function initialize!(simulator::TauLeapSimulator, state, model, tfinal)
  # reset the current system time
  simulator.t = 0.0

  simulator.tfinal = tfinal

  # set the number of jumps to zero
  fill!(simulator.next_leap_jumps, 0)

  # initialize the simulation algorithm
  initialize!(simulator.algorithm, state, model)

  # seed the stepper with the first leap
  generate_next_leap!(simulator)

  return nothing
end

# carry out the next leap
@inline function step!(simulator::TauLeapSimulator, state, model)
  # unpack information
  n = simulator.next_leap_jumps
  s = simulator.next_leap_time

  # update the system time
  simulator.t = get_new_time(simulator)

  # update the system state according to the proposed leap
  # NOTE: at this point it is assumed that the leap is valid
  simulator.execute_leap!(state, n)

  # update statistics
  # TODO

  # update the simulation algorithm given the leap event
  update!(simulator.algorithm, state, model)

  # generate the next leap
  generate_next_leap!(simulator)

  return nothing
end

@inline function generate_next_leap!(simulator::TauLeapSimulator)
  # unpack
  v = simulator.next_leap_jumps

  # set a limit on how large a leap may be
  Δt_max = simulator.tfinal - simulator.t

  # ask the algorithm for the next random leap
  v, s = generate_leap!(v, simulator.algorithm, Δt_max)

  # update the stepper's fields
  simulator.next_leap_jumps = v
  simulator.next_leap_time = s

  return nothing
end

@inline get_new_time(simulator::TauLeapSimulator) = get_new_time(simulator.algorithm, simulator.t, simulator.next_leap_time)

##### helper functions #####

@inline function no_update_step!(simulator::TauLeapSimulator, state, model)
  # unpack information
  n = simulator.next_leap_jumps
  s = simulator.next_leap_time

  # update the system time
  simulator.t = get_new_time(simulator)

  # update the system state according to the proposed leap
  # NOTE: at this point it is assumed that the leap is valid
  simulator.execute_leap!(state, n)

  return nothing
end


text:
cond1 = make_condition_dict(num_nets=750,mcts_d=[8],n_bnds=[8,120],seed=5,
                            role=:trusted)
cond2 = make_condition_dict(num_nets=250,mcts_d=[3],n_bnds=[8,120],seed=4,
                            role=:candidate)
cond3 = make_condition_dict(num_nets=250,mcts_d=[1],n_bnds=[8,120],seed=3,
                            role=:candidate)

exp_dict = Dict(:name=>"net_transition_vary", :conditions=>Dict(1=>cond1, 2=>cond2, 3=> cond3))

fldr = "logs"
for x in keys(exp_dict[:conditions])
    fname = make_useful_fname(exp_dict,x)
    exp_dict[:conditions][x][:fname] = fname
    exp_dict[:conditions][x][:fldr] = fldr
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
using WatchJuliaBurn
using Test
using LinearAlgebra

@testset "WatchJuliaBurn.jl" begin
    ## Base
    @test_throws ArgumentError c╯°□°ↄ╯(ArgumentError("Great Success"))
    @test 🗺(sin, 1:5) == sin.(1:5)
    @test 📖(:a => 2.0) == Dict(:a => 2.0)
    @test ✅
    @test 👍
    @test 👌
    @test !👎
    @test 🥩"Amazing" == raw"Amazing"
    @test 🔥(📖(:a => 2.0, :b => 3.0), :a) == 📖(:b => 3.0)
    @test_nowarn 🖨("Yes!")
    @test 📡(sin, 1:5) == sin.(1:5)
    @test ⬛ === nothing
    @test 🕵(x->x==1, [1, 2, 1]) == [1, 3]
    @test_nowarn ☝(👍)
    if VERSION >= v"1.5"
        @test ⛰(IOBuffer("Brilliant"), Char) == 'B'
    end

    ## Arrays
    @test 😻([1], [2]; dims=1) == [1, 2]
    @test ⬇😻([1], [2]) == [1, 2]
    @test ➡😻([1], [2]) == [1 2]
    @test 😹([1], [2]; dims=1) == [1, 2]
    @test ⬇😹([1], [2]) == [1, 2]
    @test ➡😹([1], [2]) == [1 2]
    @test 🐈([1], [2]; dims=1) == [1, 2]
    @test ⬇🐈([1], [2]) == [1, 2]
    @test ➡🐈([1], [2]) == [1 2]
    @test 🔢([1 0; 0 1]) == [1 0; 0 1]
    @test 🧺(1:3) == [1, 2, 3]
    if VERSION >= v"1.2"
        @eval @test_nowarn $(Symbol(Char(0x0001fa93) * Char(0x0001fa93)))(rand(3, 3)) # 🪓🪓
    end
    # @test_nowarn 🪟(rand(3, 3), 1:2, 1:2)

    ## Math
    @test 🥧 ≈ 3.1415 atol=1e-4
    @test 🍰 ≈ 3.1415 atol=1e-4
    @test 🧑🏻➡🧑🏽(2.0) == tan(2.0)
    if VERSION >= v"1.5"
        @eval @test $(Symbol(Char(0x0001fab5)))(1.0) == log(1.0) # 🪵
    end
    @test 🗜(5.0, 1.0, 2.0) == 2.0
    @test 👔💪(-2) == 2
    @test 🎽💪(-2) == 2
    @test 🛸❓(1im) == 👎
    @test 🔮(1 + 2im) == 2

    ## Monkey try/catch/finally
    include("monkeytests.jl")

end


text:
__precompile__(false)

module ReviseTest

square(x) = x^2

cube(x) = x^4  # should be x^3, but this simulates a mistake

module Internal

mult2(x) = 2*x
mult3(x) = 4*x  # oops
mult4(x) = -x

"""
This has a docstring
"""
unchanged(x) = x

unchanged2(@nospecialize(x)) = x

end  # Internal

end  # ReviseTest


text:
module MaxHelpingHandUpsideDownJumpThru

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/UpsideDownJumpThru" UpsideDownJumpThru(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, texture::String="wood",
    animationDelay::Number=0.0, pushPlayer::Bool=false, attached::Bool=false)

textures = ["wood", "dream", "temple", "templeB", "cliffside", "reflection", "core", "moon"]
const placements = Ahorn.PlacementDict(
    "Upside Down Jump Through (max480's Helping Hand)" => Ahorn.EntityPlacement(
        UpsideDownJumpThru,
        "rectangle",
        Dict{String, Any}(
            "texture" => "wood"
        )
    )
)

quads = Tuple{Integer, Integer, Integer, Integer}[
    (0, 0, 8, 7) (8, 0, 8, 7) (16, 0, 8, 7);
    (0, 8, 8, 5) (8, 8, 8, 5) (16, 8, 8, 5)
]

Ahorn.editingOptions(entity::UpsideDownJumpThru) = Dict{String, Any}(
    "texture" => textures
)

Ahorn.minimumSize(entity::UpsideDownJumpThru) = 8, 0
Ahorn.resizable(entity::UpsideDownJumpThru) = true, false

function Ahorn.selection(entity::UpsideDownJumpThru)
    x, y = Ahorn.position(entity)
    width = Int(get(entity.data, "width", 8))

    return Ahorn.Rectangle(x, y, width, 8)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::UpsideDownJumpThru, room::Maple.Room)
    texture = get(entity.data, "texture", "wood")
    texture = texture == "default" ? "wood" : texture

    # Values need to be system specific integer
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 8))

    startX = div(x, 8) + 1
    stopX = startX + div(width, 8) - 1
    startY = div(y, 8) + 1
    animated = Number(get(entity.data, "animationDelay", 0)) > 0

    Ahorn.Cairo.save(ctx)

    Ahorn.scale(ctx, 1, -1)

    len = stopX - startX
    for i in 0:len
        if animated
            Ahorn.drawImage(ctx, "objects/jumpthru/$(texture)00", 8 * i, -8)
        else
            connected = false
            qx = 2
            if i == 0
                connected = get(room.fgTiles.data, (startY, startX - 1), false) != '0'
                qx = 1

            elseif i == len
                connected = get(room.fgTiles.data, (startY, stopX + 1), false) != '0'
                qx = 3
            end

            quad = quads[2 - connected, qx]
            Ahorn.drawImage(ctx, "objects/jumpthru/$(texture)", 8 * i, -8, quad...)
        end
    end

    Ahorn.Cairo.restore(ctx)
end

end

max_stars_repo_name: swoolcock/MaxHelpingHand

max_issues_repo_name: swoolcock/MaxHelpingHand

max_forks_repo_name: swoolcock/MaxHelpingHand


text:
julia> A = [2 9 4; 7 5 3; 6 1 8]
julia> det(A), perm(A)
(-360.0,900)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# This file was generated, do not modify it. # hide
ŷ = exp(ẑ);


text:
cd(@__DIR__)
using AdvancedPS, Random, KernelDensity, Distributions
include("KDE.jl")

function loglike(μ, σ, data)
    simdata = rand(Normal(μ, σ), 10_000)
    kd = kernel(simdata)
    dist = InterpKDE(kd)
    #println("proposal: ",μ," ",σ)
    like = max.(1e-10, pdf(dist, data))
    return sum(log.(like))
end

priors = (
    μ=(Normal(0, 10),),
    σ=(Truncated(Cauchy(0, 1), 0.0, Inf),)
)

bounds = ((-Inf,Inf),(0.0,Inf))

data = rand(Normal(0.0, 1.0), 50)

model = DEModel(priors=priors, model=x->loglike(x..., data))

de = DE(bounds=bounds, visualize=false, burnin=1000, priors=priors, progress=true)
n_iter = 2000
chains = sample(model, de, n_iter)
println(chains)

max_stars_repo_path: DEMCMC_Examples/KDE Example.jl

max_issues_repo_path: DEMCMC_Examples/KDE Example.jl

max_forks_repo_path: DEMCMC_Examples/KDE Example.jl


text:
function forward_pass!(p_data::PolicyData, m_data::ModelData, s_data::SolverData;
    max_iter = 25)

    # reset solver status
    s_data.status = false

    # previous cost
    J_prev = s_data.obj

    # gradient of Lagrangian
    lagrangian_gradient!(s_data, p_data, m_data)

    # line search with rollout
    s_data.α = 1.0
    iter = 1
    while true
        iter > max_iter && (@error "forward pass failure", break)

        J = Inf
        try
            rollout!(p_data, m_data, α = s_data.α)
            J = objective!(s_data, m_data, mode = :current)
            Δz!(m_data)
        catch
            @warn "rollout failure"
            fill!(m_data.z, 0.0)
        end

        if J < J_prev + 0.001 * s_data.α * s_data.gradient' * m_data.z
            # update nominal
            m_data.x̄ .= deepcopy(m_data.x)
            m_data.ū .= deepcopy(m_data.u)
            s_data.obj = J
            s_data.status = true
            break
        else
            s_data.α *= 0.5
            iter += 1
        end
    end
end


text:
include("Graph.jl")
include("Centrality.jl")

fileName = string("Dataset/", ARGS[1], ".txt")
G0 = readGraph(fileName, "unweighted")
G = getLLC(G0)
lg = open("log.txt", "a")
println(lg, ARGS[1], " ", G0.n, " ", G0.m, " ", G.n, " ", G.m)

println(lg, "DC : ", getPc(Degree(G)))
println(lg, "BC : ", getPc(BetweennessCentrality(G)))
println(lg, "CC : ", getPc(ClosenessCentrality(G)))
println(lg, "EC : ", getPc(EigenvectorCentrality(G)))
println(lg, "PR : ", getPc(PageRank(G)))
println(lg, "IC : ", getPc(InformationCentrality(G)))
println(lg, "FDC : ", getPc(ForestDistanceClosenessCentrality(G)))
println(lg)


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ e6db292c-d2be-11eb-1c25-05ed50a9256c
using FourierTools, TestImages, ImageShow, Plots, Colors, FFTW, Noise, PlutoUI

# ╔═╡ 4958079f-da85-4357-94df-b122cf01c958
FFTW.set_num_threads(12)

# ╔═╡ 2a47b7a9-7767-44a7-8dda-480ba1b41f57
md"$\mathcal{F}[U](\mathbf k) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} U(\mathbf r) \cdot \exp[-i  (\mathbf x \cdot \mathbf k)] \, \mathrm d \mathbf x$

$\mathcal{F}^{-1}[\tilde U](\mathbf x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} U(\mathbf r) \cdot \exp[i  (\mathbf x \cdot \mathbf k)] \, \mathrm d \mathbf k$
";

# ╔═╡ b71a5e39-0805-4e3f-aa83-c3972f6c54af
begin
	img = Float32.(testimage("fabio_gray"));
	img_1D = img[:, 200]
end;

# ╔═╡ c39440fc-17e6-440c-9e10-7b2c9a7e9331
md"### FourierTools.jl | Working with the Frequency Space

* **Felix Wechsler**: Leibniz Institute of Photonic Technology in Jena, Germany
* **Rainer Heintzmann**: Leibniz Institute of Photonic Technology in Jena, Germany
"

# ╔═╡ 6e445cd5-6e51-4a8f-ad5c-2636b3da68f4
md"## FourierTools.jl

* A package built on FFTW.jl
* Provides several helper functions
* implements algorithms (like interpolation) based on FFTs
"

# ╔═╡ e28dc936-04da-48cb-b324-73632d008290
md"
### Definition of Fourier Transform

Complex valued sum

$\mathcal F[x_k] = X_k = \sum_{n=0}^{N-1} x_n \cdot \exp\left[-\frac{i 2 \pi}{N} k n\right] =$

$=\sum_{n=0}^{N-1} x_n \cdot \left[ \cos\left(\frac{2 \pi}{N} k n\right) - i \cdot \sin\left(\frac{2 \pi}{N} k n\right)\right]$

Inverse transform

$\mathcal F^{-1}[X_k] = x_k = \frac{1}{N} \sum_{k=0}^{N-1} X_k \cdot \exp\left[\frac{i 2 \pi}{N} k n\right]$

* it decomposes a signal into harmonic basis functions
* the signal is represented in the _Fourier_ or _Frequency_ space

### Fast Fourier Transform - FFT
* naive sum calculation takes $\mathcal O(N^2)$ operations for a 1D dataset with $N$ points
* FFT algorithms evaluates this with $\mathcal O(N \cdot \log(N))$ operations
* In Julia available _FFTW.jl_
"

# ╔═╡ e8f1e951-26dd-4c0e-83f9-4a2c60b0ebed
md" ## FFT of an dataset with FFTW.jl"

# ╔═╡ 3adab086-0c77-40e5-990a-6b89b183cc4d
begin
	x = range(0, 50, length=999)
	f(x) = (0.1 * sin(17.3 + x) + 0.5 * cos(1.2*x) + 0.3 * sin(5*x)^2 + 0.05 * sin(13.1*x) + sin(27.2 * x)) * exp(-x/20)
	data = add_gauss(f.(x), 0.1)
end;

# ╔═╡ 850869e6-32da-48cf-b2a0-081cc9f8e524
plot(x, data, title="function to transform")

# ╔═╡ d6f9bd7a-ddcb-40b5-8509-f4dac3f8bd5d
begin
	plot(real(fft(data)), title="FFT of the function", xlabel="frequency", label="real part")
	plot!(imag(fft(data)), label="imaginary part")
end

# ╔═╡ b4b57f32-e405-4463-b0fa-625e9d54df34
md"
* Data: $[x_1, x_2, ..., x_N]$
* FFT of Data: $[f_0, f_1, ..., f_{N/2}, f_{-N/2}, f_{-N/2+1}, ...,f_{-1}]$
* (for real data negative and positive frequencies are the complex conjugate)
"

# ╔═╡ f0746bc6-f50a-49f1-96d2-0379498a0168
begin
	plot(fftshift(fftfreq(size(data, 1))), real(ffts(data)), xlabel="frequency", label="real part", title="FFT of the dataset but shifted")
	plot!(fftshift(fftfreq(size(data, 1))), imag(ffts(data)), label="imaginary part")
end

# ╔═╡ 979d641e-83ad-4484-bc28-912947711936
fftshift(fft(data)) == ffts(data)

# ╔═╡ 0ba41c40-21a3-4e26-99fe-cacd76d47a67
typeof(ffts(data))

# ╔═╡ cc160200-5911-4c80-b651-a54e1608f186
begin
	N_low = 64
	x_min = 0.0
	x_max = 8*2π

	xs_low = range(x_min, x_max, length=N_low+1)[1:N_low]
	xs_high = range(x_min, x_max, length=5000)[1:end-1]
	f2(x) = sin(0.5*x) + cos(x) + cos(2 * x) + sin(0.25*x)
	arr_low = f2.(xs_low)
	arr_high = f2.(xs_high)
end;

# ╔═╡ 7acf4c37-397d-4c1f-8735-022888390cbe
md"## Different FFT methods
* `fft` interpretes the data as top-left centered

"

# ╔═╡ b15df267-6606-4583-8ce9-3dddf1ecf3f9
fft([1.0,0,0,0,0])

# ╔═╡ 76bb0c74-1a52-47ff-adf2-862cfea15da5
fft([0,0, 1.0 ,0,0]) # single delta peak

# ╔═╡ ec839c01-a62f-46eb-8e25-15d04a07a79f
ft([0,0, 1.0 ,0,0])

# which is equivalent to fftshift(fft(ifftshift([0,0, 1.0 ,0,0])))

# ╔═╡ c7f0e20e-7c80-4a47-abdf-800f26ca999a
ffts([0,0, 1.0 ,0,0])

# which is equivalent to fftshift(fft([0,0, 1.0 ,0,0]))

# ╔═╡ a725ae98-e2cc-439c-9cc1-692972ac98e7
md"## Calculate sinc interpolation

* Fourier coefficients represent coefficients of Fourier series
* Inverse transform is essentially a re-evaluation of a Fourier Series
* Evaluating this series at more samples points results in an interpolated signal
"

# ╔═╡ d5717c9c-2c16-47da-8533-a28acf189e55
begin
	N = 1000
	xs_interp = range(x_min, x_max, length=N+1)[1:N]
	arr_interp = resample(arr_low, N)
end;

# ╔═╡ 3d8f11a3-4ebf-4d36-a227-b2cafa13e1d6
begin
	scatter(xs_low, arr_low, legend=:bottomleft, markersize=2, label="Low sampling")
	plot!(xs_interp, arr_interp, label="FFT based sinc interpolation", linestyle=:dash)
	plot!(xs_high, arr_high, linestyle=:dashdotdot, label="High sampling")
end

# ╔═╡ 551f179b-46ad-425f-9dd6-d8f6469a09cf
[Gray.(img) Gray.(resample(img, (256, 712)))]

# ╔═╡ a5977c5c-dbcd-402d-b7a3-56614968ccb7
begin
	N_ds = 32
	xs_ds = range(x_min, x_max, length=N_ds+1)[1:N_ds]
	arr_ds = real(resample(arr_high, (N_ds,)))
end;

# ╔═╡ 7f3be670-cfe3-4e60-b38f-ebacb4eaf8d3
md"## Fourier Shift Theorem

$\mathcal{F}[U(x+ Δx)](\mathbf k) = \mathcal{F}[U(x)](\mathbf k) \exp(i k \Delta x)$

$\Longrightarrow$
$U(x+ Δx) = \mathcal{F}^{-1}\bigg[\mathcal{F}[U(x)](\mathbf k) \exp(i k \Delta x)\bigg](x)$

* which allows for efficient sub-pixel shifting
"

# ╔═╡ 1ab61f99-799f-4c34-b0d1-205c01c20dbc
@bind offset Slider(-10:0.01:10, show_value=true)

# ╔═╡ 55667950-4542-4505-9d21-e44ea8b6074d
begin
	f3(x) = cos(4π * x / 30)
	x1 = 1:30
	x2 = x1 .+ 3
	y1 = f3.(x1)
	y2 = f3.(x2)
	y3 = shift(y2, tuple(offset))
end;

# ╔═╡ 8fc68dfc-0fb7-4a5a-bef1-ae0b33ddf4f6
shift(y2, (offset,))

# ╔═╡ cfb02910-fb2d-4de0-bc81-8ecea50e5cb5
begin
	plot(y1, label="Original signal")
	plot!(y2, label="Shifted signal")
	plot!(y3, label="Fourier shifted with $offset")
end

# ╔═╡ 105344de-8617-47fd-aeed-84102f19eaeb
md"## Based on Shifting $\Longrightarrow$ Rotation

"

# ╔═╡ 344f9221-4213-4540-bb40-8a2de2d9b48c
begin
	z = zeros(Float32, round.(Ref(Int), 1.5 .* size(img)))
	FourierTools.center_set!(z, img)
end;

# ╔═╡ 27a7bcae-e50f-4225-bd7f-45934ec55be3
@bind ϕ Slider(-180:1:180, show_value=true)

# ╔═╡ ac331005-48ca-4be8-b0a8-85f6f740f09a
FourierTools.rotate(z, ϕ);

# ╔═╡ bf2c54f1-8110-4e25-b89c-7276cf55f42a
Gray.(FourierTools.rotate(z, ϕ))

# ╔═╡ a5ea8226-ff14-432c-974e-68e1a79155cf
Gray.(FourierTools.rotate(img, ϕ))

# ╔═╡ e1600b0a-b8ef-49b8-8ed0-f9c4ae0c436b

# ╔═╡ 7ff8de2d-85fd-402a-88eb-16427e1e0b47
md"## Convolution
Efficient FFT convolutions for both real and complex functions are included as well
"

# ╔═╡ 03effebd-e918-4847-b13a-7d41d4f159e7
md"Blurring strength"

# ╔═╡ e9862df6-8d1a-4c5e-88ab-1da28c2652af
@bind pow Slider(1:0.01:6, show_value=true)

# ╔═╡ bf631ded-b54c-4768-9a8d-825ab7612695
begin
	x3 = fftpos(1, size(img, 1))
	kernel = exp.(.- 10.0 .^pow .* (x3 .^2 .+ x3'.^2))
	kernel ./= sum(kernel)
end;

# ╔═╡ f4870819-9db6-4bdc-8c8d-141c8fa67461
img_blurry = conv_psf(img, kernel);

# ╔═╡ 9d22e799-fe80-47ef-a326-1f5f1ccc07f8
[Gray.(img) Gray.(img_blurry)]

# ╔═╡ c530d7ea-60bd-4458-b3ca-80b509fb4f33
md"# Conclusion

* FourierTools.jl provides convenient wrappers to the FFTW library
* offers efficient algorithms like sinc interpolation based on FFTs
* deals with the intricacies of the correct handling of even/odd sized arrays in Fourier space
"

# ╔═╡ Cell order:
# ╠═e6db292c-d2be-11eb-1c25-05ed50a9256c
# ╠═4958079f-da85-4357-94df-b122cf01c958
# ╟─2a47b7a9-7767-44a7-8dda-480ba1b41f57
# ╟─b71a5e39-0805-4e3f-aa83-c3972f6c54af
# ╟─c39440fc-17e6-440c-9e10-7b2c9a7e9331
# ╟─6e445cd5-6e51-4a8f-ad5c-2636b3da68f4
# ╟─e28dc936-04da-48cb-b324-73632d008290
# ╟─e8f1e951-26dd-4c0e-83f9-4a2c60b0ebed
# ╟─3adab086-0c77-40e5-990a-6b89b183cc4d
# ╟─850869e6-32da-48cf-b2a0-081cc9f8e524
# ╟─d6f9bd7a-ddcb-40b5-8509-f4dac3f8bd5d
# ╟─b4b57f32-e405-4463-b0fa-625e9d54df34
# ╟─f0746bc6-f50a-49f1-96d2-0379498a0168
# ╠═979d641e-83ad-4484-bc28-912947711936
# ╠═0ba41c40-21a3-4e26-99fe-cacd76d47a67
# ╟─cc160200-5911-4c80-b651-a54e1608f186
# ╟─7acf4c37-397d-4c1f-8735-022888390cbe
# ╠═b15df267-6606-4583-8ce9-3dddf1ecf3f9
# ╠═76bb0c74-1a52-47ff-adf2-862cfea15da5
# ╠═ec839c01-a62f-46eb-8e25-15d04a07a79f
# ╠═c7f0e20e-7c80-4a47-abdf-800f26ca999a
# ╟─a725ae98-e2cc-439c-9cc1-692972ac98e7
# ╟─d5717c9c-2c16-47da-8533-a28acf189e55
# ╟─3d8f11a3-4ebf-4d36-a227-b2cafa13e1d6
# ╠═551f179b-46ad-425f-9dd6-d8f6469a09cf
# ╟─a5977c5c-dbcd-402d-b7a3-56614968ccb7
# ╟─7f3be670-cfe3-4e60-b38f-ebacb4eaf8d3
# ╟─1ab61f99-799f-4c34-b0d1-205c01c20dbc
# ╠═8fc68dfc-0fb7-4a5a-bef1-ae0b33ddf4f6
# ╟─55667950-4542-4505-9d21-e44ea8b6074d
# ╟─cfb02910-fb2d-4de0-bc81-8ecea50e5cb5
# ╟─105344de-8617-47fd-aeed-84102f19eaeb
# ╟─344f9221-4213-4540-bb40-8a2de2d9b48c
# ╟─27a7bcae-e50f-4225-bd7f-45934ec55be3
# ╠═ac331005-48ca-4be8-b0a8-85f6f740f09a
# ╟─bf2c54f1-8110-4e25-b89c-7276cf55f42a
# ╟─a5ea8226-ff14-432c-974e-68e1a79155cf
# ╠═e1600b0a-b8ef-49b8-8ed0-f9c4ae0c436b
# ╟─7ff8de2d-85fd-402a-88eb-16427e1e0b47
# ╟─03effebd-e918-4847-b13a-7d41d4f159e7
# ╟─e9862df6-8d1a-4c5e-88ab-1da28c2652af
# ╠═f4870819-9db6-4bdc-8c8d-141c8fa67461
# ╠═9d22e799-fe80-47ef-a326-1f5f1ccc07f8
# ╟─bf631ded-b54c-4768-9a8d-825ab7612695
# ╟─c530d7ea-60bd-4458-b3ca-80b509fb4f33


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SCIP_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SCIP")
JLLWrappers.@generate_main_file("SCIP", UUID("e5ac4fe4-a920-5659-9bf8-f9f73e9e79ce"))
end  # module SCIP_jll


text:
cachetypeof(::Type{<:CentralizedRecycler}) = CentralizedCache

# Workaround for avoid storing types as `DataType``
CentralizedRecycler(::Type{Factory}, cache::CentralizedCache{T}) where {Factory,T} =
    CentralizedRecycler{T,Type{Factory}}(Factory, cache)
CentralizedRecycler{T}(::Type{Factory}, cache::CentralizedCache{T}) where {Factory,T} =
    CentralizedRecycler{T,Type{Factory}}(Factory, cache)

CentralizedRecycler{T}(factory::Factory, cache::CentralizedCache{T}) where {Factory,T} =
    CentralizedRecycler{T,Factory}(factory, cache)

CentralizedCache{T}() where {T} = _CentralizedCache(T, nothing)

function CentralizedCache{T}(objects) where {T}
    node = nothing
    for x in objects
        node = StackNode{T}(x, node)
    end
    return _CentralizedCache(T, node)
end

function Recyclers.maybeget!(cache::CentralizedCache{T}) where {T}
    top = @atomic :acquire cache.top
    while true
        old = if top === nothing
            return nothing
        else
            top
        end::StackNode{T}
        next = old.next
        (top, success) = @atomicreplace :acquire_release :acquire cache.top old => next
        success && return Some{T}(old.object)
    end
end

function Recyclers.recycle!(cache::CentralizedCache{T}, object::T) where {T}
    new = StackNode{T}(object, nothing)
    top = @atomic :monotonic cache.top
    while true
        new.next = top
        (top, success) = @atomicreplace :acquire_release :acquire cache.top top => new
        success && return true
    end
end

function Base.empty!(cache::CentralizedCache)
    @atomic :release cache.top = nothing
    return cache
end

function Recyclers.unsafe_empty!(cache::CentralizedCache)
    @atomic :monotonic cache.top = nothing
    return cache
end

Recyclers.unsafe_destroy!(cache::CentralizedCache) = Recyclers.unsafe_empty!(cache)
Recyclers.unsafe_init!(cache::CentralizedCache) = Recyclers.unsafe_empty!(cache)


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

#proj
ilab = 1:10
olab = [randstring(1) for i =1:10]
fst = WFST(ilab,olab)
add_arc!(fst, 1, 1, ilab[1] , olab[1] , rand())
add_arc!(fst, 1, 2, ilab[2] , olab[2] , rand())
add_arc!(fst, 2, 3, ilab[3] , olab[3] , rand())
add_arc!(fst, 2, 2, ilab[4] , olab[4] , rand())
add_arc!(fst, 3, 3, ilab[5] , olab[5] , rand())
add_arc!(fst, 3, 4, ilab[6] , olab[6] , rand())
add_arc!(fst, 4, 4, ilab[7] , olab[7] , rand())
add_arc!(fst, 4, 5, ilab[8] , olab[8] , rand())
add_arc!(fst, 5, 5, ilab[9] , olab[9] , rand())
add_arc!(fst, 5, 6, ilab[10], olab[10], rand())
initial!(fst,1)
final!(fst,6,rand())

pfst = proj(get_ilabel,fst)
@test isfinal(pfst,6)
@test isinitial(pfst,1)
@test is_acceptor(pfst)
@test get_weight(fst,6) == get_weight(pfst,6)
@test all([a in ilab for a in get_oalphabet(pfst)])
@test all([a in ilab for a in get_ialphabet(pfst)])

pfst = proj(get_olabel,fst)
@test isfinal(pfst,6)
@test isinitial(pfst,1)
@test is_acceptor(pfst)
@test get_weight(fst,6) == get_weight(pfst,6)
@test all([a in olab for a in get_oalphabet(pfst)])
@test all([a in olab for a in get_ialphabet(pfst)])


text:
using LightGraphs
using MetaGraphs

"""
    disconnect_param!(obj::AbstractCompositeComponentDef, comp_def::AbstractComponentDef, param_name::Symbol)

Remove any parameter connections for a given parameter `param_name` in a given component
`comp_def` which must be a direct subcomponent of composite `obj`.
"""
function disconnect_param!(obj::AbstractCompositeComponentDef, comp_def::AbstractComponentDef, param_name::Symbol)
    # If the path isn't set yet, we look for a comp in the eventual location
    path = @or(comp_def.comp_path, ComponentPath(obj, comp_def.name))

    # @info "disconnect_param!($(obj.comp_path), $path, :$param_name)"

    if is_descendant(obj, comp_def) === nothing
        error("Cannot disconnect a component ($path) that is not within the given composite ($(obj.comp_path))")
    end

    filter!(x -> !(x.dst_comp_path == path && x.dst_par_name == param_name), obj.internal_param_conns)

    if obj isa ModelDef
        filter!(x -> !(x.comp_path == path && x.param_name == param_name), obj.external_param_conns)
    end
    dirty!(obj)
end

"""
    disconnect_param!(obj::AbstractCompositeComponentDef, comp_name::Symbol, param_name::Symbol)

Remove any parameter connections for a given parameter `param_name` in a given component
`comp_def` which must be a direct subcomponent of composite `obj`.
"""
function disconnect_param!(obj::AbstractCompositeComponentDef, comp_name::Symbol, param_name::Symbol)
    comp = compdef(obj, comp_name)
    comp === nothing && error("Did not find $comp_name in composite $(printable(obj.comp_path))")
    disconnect_param!(obj, comp, param_name)
end

"""
    disconnect_param!(obj::AbstractCompositeComponentDef, comp_path::ComponentPath, param_name::Symbol)

Remove any parameter connections for a given parameter `param_name` in the component identified by
`comp_path` which must be under the composite `obj`.
"""
function disconnect_param!(obj::AbstractCompositeComponentDef, comp_path::ComponentPath, param_name::Symbol)
    if (comp_def = find_comp(obj, comp_path)) === nothing
        return
    end
    disconnect_param!(obj, comp_def, param_name)
end

# Default string, string unit check function
verify_units(unit1::AbstractString, unit2::AbstractString) = (unit1 == unit2)

function _check_labels(obj::AbstractCompositeComponentDef,
                       comp_def::AbstractComponentDef, param_name::Symbol, ext_param::ArrayModelParameter)
    param_def = parameter(comp_def, param_name)

    t1 = eltype(ext_param.values)
    t2 = eltype(param_def.datatype)
    if !(t1 <: Union{Missing, t2})
        error("Mismatched datatype of parameter connection: Component: $(comp_def.comp_id) ($t1), Parameter: $param_name ($t2)")
    end

    comp_dims  = dim_names(param_def)
    param_dims = dim_names(ext_param)

    if ! isempty(param_dims) && size(param_dims) != size(comp_dims)
        d1 = size(comp_dims)
        d2 = size(param_dims)
        error("Mismatched dimensions of parameter connection: Component: $(comp_def.comp_id) ($d1), Parameter: $param_name ($d2)")
    end

    # Don't check sizes for ConnectorComps since they won't match.
    if nameof(comp_def) in (:ConnectorCompVector, :ConnectorCompMatrix)
        return nothing
    end

    # index_values = indexvalues(obj)

    for (i, dim) in enumerate(comp_dims)
        if isa(dim, Symbol)
            param_length = size(ext_param.values)[i]
            if dim == :time
                t = dimension(obj, :time)
                first = find_first_period(comp_def)
                last = find_last_period(comp_def)
                comp_length = t[last] - t[first] + 1
            else
                comp_length = dim_count(obj, dim)
            end
            if param_length != comp_length
                error("Mismatched data size for a parameter connection: dimension :$dim in $(comp_def.comp_id) has $comp_length elements; external parameter :$param_name has $param_length elements.")
            end
        end
    end
end

"""
    connect_param!(obj::AbstractCompositeComponentDef, comp_name::Symbol, param_name::Symbol, ext_param_name::Symbol;
                   check_labels::Bool=true)

Connect a parameter `param_name` in the component `comp_name` of composite `obj` to
the external parameter `ext_param_name`.
"""
function connect_param!(obj::AbstractCompositeComponentDef, comp_name::Symbol,
                        param_name::Symbol, ext_param_name::Symbol;
                        check_labels::Bool=true)
    comp_def = compdef(obj, comp_name)
    connect_param!(obj, comp_def, param_name, ext_param_name, check_labels=check_labels)
end

function connect_param!(obj::AbstractCompositeComponentDef, comp_def::AbstractComponentDef,
                        param_name::Symbol, ext_param_name::Symbol; check_labels::Bool=true)
    ext_param = external_param(obj, ext_param_name)

    if ext_param isa ArrayModelParameter && check_labels
        _check_labels(obj, comp_def, param_name, ext_param)
    end

    disconnect_param!(obj, comp_def, param_name)    # calls dirty!()

    comp_path = @or(comp_def.comp_path, ComponentPath(obj.comp_path, comp_def.name))
    conn = ExternalParameterConnection(comp_path, param_name, ext_param_name)
    add_external_param_conn!(obj, conn)

    return nothing
end

"""
    connect_param!(obj::AbstractCompositeComponentDef,
        dst_comp_path::ComponentPath, dst_par_name::Symbol,
        src_comp_path::ComponentPath, src_var_name::Symbol,
        backup::Union{Nothing, Array}=nothing;
        ignoreunits::Bool=false, offset::Int=0)

Bind the parameter `dst_par_name` of one component `dst_comp_path` of composite `obj` to a
variable `src_var_name` in another component `src_comp_path` of the same model using
`backup` to provide default values and the `ignoreunits` flag to indicate the need to
check match units between the two.  The `offset` argument indicates the offset between
the destination and the source ie. the value would be `1` if the destination component
parameter should only be calculated for the second timestep and beyond.
"""
function connect_param!(obj::AbstractCompositeComponentDef,
                        dst_comp_path::ComponentPath, dst_par_name::Symbol,
                        src_comp_path::ComponentPath, src_var_name::Symbol,
                        backup::Union{Nothing, Array}=nothing;
                        ignoreunits::Bool=false, offset::Int=0)

    # remove any existing connections for this dst parameter
    disconnect_param!(obj, dst_comp_path, dst_par_name)  # calls dirty!()

    dst_comp_def = compdef(obj, dst_comp_path)
    src_comp_def = compdef(obj, src_comp_path)

    # @info "dst_comp_def: $dst_comp_def"
    # @info "src_comp_def: $src_comp_def"

    if backup !== nothing
        # If value is a NamedArray, we can check if the labels match
        if isa(backup, NamedArray)
            dims = dimnames(backup)
            check_parameter_dimensions(obj, backup, dims, dst_par_name)
        else
            dims = nothing
        end

        # Check that the backup data is the right size
        if size(backup) != datum_size(obj, dst_comp_def, dst_par_name)
            error("Cannot connect parameter; the provided backup data is the wrong size. ",
                  "Expected size $(datum_size(obj, dst_comp_def, dst_par_name)) but got $(size(backup)).")
        end

        # some other check for second dimension??
        dst_param = parameter(dst_comp_def, dst_par_name)
        dst_dims  = dim_names(dst_param)

        # convert number type and, if it's a NamedArray, convert to Array
        backup = convert(Array{Union{Missing, number_type(obj)}}, backup)
        first = first_period(obj, dst_comp_def)

        T = eltype(backup)

        dim_count = length(dst_dims)

        if dim_count == 0
            values = backup
        else
            ti = get_time_index_position(dst_param)

            if isuniform(obj)
                # use the first from the comp_def not the ModelDef
                stepsize = step_size(obj)
                last = last_period(obj, dst_comp_def)
                values = TimestepArray{FixedTimestep{first, stepsize, last}, T, dim_count, ti}(backup)
            else
                times = time_labels(obj)
                # use the first from the comp_def
                first_index = findfirst(isequal(first), times)
                values = TimestepArray{VariableTimestep{(times[first_index:end]...,)}, T, dim_count, ti}(backup)
            end

        end

        set_external_array_param!(obj, dst_par_name, values, dst_dims)
        backup_param_name = dst_par_name

    else
        # If backup not provided, make sure the source component covers the span of the destination component
        src_first, src_last = first_and_last(src_comp_def)
        dst_first, dst_last = first_and_last(dst_comp_def)
        if (dst_first !== nothing && src_first !== nothing && dst_first < src_first) ||
            (dst_last  !== nothing && src_last  !== nothing && dst_last  > src_last)
            src_first = printable(src_first)
            src_last  = printable(src_last)
            dst_first = printable(dst_first)
            dst_last  = printable(dst_last)
            error("""Cannot connect parameter: $src_comp_path runs only from $src_first to $src_last,
whereas $dst_comp_path runs from $dst_first to $dst_last. Backup data must be provided for missing years.
Try calling:
    `connect_param!(m, comp_name, par_name, comp_name, var_name, backup_data)`""")
        end

        backup_param_name = nothing
    end

    # Check the units, if provided
    if ! ignoreunits && ! verify_units(variable_unit(src_comp_def, src_var_name),
                                       parameter_unit(dst_comp_def, dst_par_name))
        error("Units of $src_comp_path:$src_var_name do not match $dst_comp_path:$dst_par_name.")
    end

    conn = InternalParameterConnection(src_comp_path, src_var_name, dst_comp_path, dst_par_name,
                                       ignoreunits, backup_param_name, offset=offset)
    add_internal_param_conn!(obj, conn)

    return nothing
end

function connect_param!(obj::AbstractCompositeComponentDef,
                        dst_comp_name::Symbol, dst_par_name::Symbol,
                        src_comp_name::Symbol, src_var_name::Symbol,
                        backup::Union{Nothing, Array}=nothing; ignoreunits::Bool=false, offset::Int=0)
    connect_param!(obj, ComponentPath(obj, dst_comp_name), dst_par_name,
                        ComponentPath(obj, src_comp_name), src_var_name,
                        backup; ignoreunits=ignoreunits, offset=offset)
end

"""
    connect_param!(obj::AbstractCompositeComponentDef,
        dst::Pair{Symbol, Symbol}, src::Pair{Symbol, Symbol},
        backup::Union{Nothing, Array}=nothing;
        ignoreunits::Bool=false, offset::Int=0)

Bind the parameter `dst[2]` of one component `dst[1]` of composite `obj`
to a variable `src[2]` in another component `src[1]` of the same composite
using `backup` to provide default values and the `ignoreunits` flag to indicate the need
to check match units between the two.  The `offset` argument indicates the offset
between the destination and the source ie. the value would be `1` if the destination
component parameter should only be calculated for the second timestep and beyond.
"""
function connect_param!(obj::AbstractCompositeComponentDef,
                        dst::Pair{Symbol, Symbol}, src::Pair{Symbol, Symbol},
                        backup::Union{Nothing, Array}=nothing; ignoreunits::Bool=false, offset::Int=0)
    connect_param!(obj, dst[1], dst[2], src[1], src[2], backup; ignoreunits=ignoreunits, offset=offset)
end

"""
    split_datum_path(obj::AbstractCompositeComponentDef, s::AbstractString)

Split a string of the form "/path/to/component:datum_name" into the component path,
`ComponentPath(:path, :to, :component)` and name `:datum_name`.
"""
function split_datum_path(obj::AbstractCompositeComponentDef, s::AbstractString)
    elts = split(s, ":")
    length(elts) != 2 && error("Can't split datum path '$s' into ComponentPath and datum name")
    return (ComponentPath(obj, elts[1]), Symbol(elts[2]))
end

"""
    connection_refs(obj::AbstractCompositeComponentDef)

Return a vector of UnnamedReference's to parameters from subcomponents that are either found in
internal connections or that have been already imported in a parameter definition.
"""
function connection_refs(obj::AbstractCompositeComponentDef)
    refs = UnnamedReference[]

    for conn in obj.internal_param_conns
        push!(refs, UnnamedReference(conn.dst_comp_path.names[end], conn.dst_par_name))
    end

    for item in values(obj.namespace)
        if item isa CompositeParameterDef
            for ref in item.refs
                push!(refs, ref)
            end
        end
    end

    return refs
end

"""
    connection_refs(obj::ModelDef)

Return a vector of UnnamedReference's to parameters from subcomponents that are either found in
internal connections or that have been already connected to external parameter values.
"""
function connection_refs(obj::ModelDef)
    refs = UnnamedReference[]

    for conn in obj.internal_param_conns
        push!(refs, UnnamedReference(conn.dst_comp_path.names[end], conn.dst_par_name))
    end

    for conn in obj.external_param_conns
        push!(refs, UnnamedReference(conn.comp_path.names[end], conn.param_name))
    end

    return refs
end

"""
    unconnected_params(obj::AbstractCompositeComponentDef)

Return a list of UnnamedReference's to parameters that have not been connected
to a value.
"""
function unconnected_params(obj::AbstractCompositeComponentDef)
    return setdiff(subcomp_params(obj), connection_refs(obj))
end

"""
    set_leftover_params!(m::Model, parameters::Dict)

Set all of the parameters in model `m` that don't have a value and are not connected
to some other component to a value from a dictionary `parameters`. This method assumes
the dictionary keys are strings that match the names of unset parameters in the model.
"""
function set_leftover_params!(md::ModelDef, parameters::Dict{T, Any}) where T
    for param_ref in unconnected_params(md)
        param_name = param_ref.datum_name
        comp_name = param_ref.comp_name
        comp_def = find_comp(md, comp_name)
        param_def = comp_def[param_name]

        # Only set the unconnected parameter if it doesn't have a default
        if param_def.default === nothing
            # check whether we need to create the external parameter
            if external_param(md, param_name, missing_ok=true) === nothing
                if haskey(parameters, string(param_name))
                    value = parameters[string(param_name)]
                    param_dims = parameter_dimensions(md, comp_name, param_name)

                    set_external_param!(md, param_name, value; param_dims = param_dims)
                else
                    error("Cannot set parameter :$param_name, not found in provided dictionary and no default value detected.")
                end
            end
            connect_param!(md, comp_name, param_name, param_name)
        end
    end
    nothing
end

# Find internal param conns to a given destination component
function internal_param_conns(obj::AbstractCompositeComponentDef, dst_comp_path::ComponentPath)
    return filter(x->x.dst_comp_path == dst_comp_path, internal_param_conns(obj))
end

function internal_param_conns(obj::AbstractCompositeComponentDef, comp_name::Symbol)
    return internal_param_conns(obj, ComponentPath(obj.comp_path, comp_name))
end

function add_internal_param_conn!(obj::AbstractCompositeComponentDef, conn::InternalParameterConnection)
    push!(obj.internal_param_conns, conn)
    dirty!(obj)
end

#
# These should all take ModelDef instead of AbstractCompositeComponentDef as 1st argument
#

# Find external param conns for a given comp
function external_param_conns(obj::ModelDef, comp_path::ComponentPath)
    return filter(x -> x.comp_path == comp_path, external_param_conns(obj))
end

function external_param_conns(obj::ModelDef, comp_name::Symbol)
    return external_param_conns(obj, ComponentPath(obj.comp_path, comp_name))
end

function external_param(obj::ModelDef, name::Symbol; missing_ok=false)
    haskey(obj.external_params, name) && return obj.external_params[name]

    missing_ok && return nothing

    error("$name not found in external parameter list")
end

function add_external_param_conn!(obj::ModelDef, conn::ExternalParameterConnection)
    push!(obj.external_param_conns, conn)
    dirty!(obj)
end

function set_external_param!(obj::ModelDef, name::Symbol, value::ModelParameter)
    # if haskey(obj.external_params, name)
    #     @warn "Redefining external param :$name in $(obj.comp_path) from $(obj.external_params[name]) to $value"
    # end
    obj.external_params[name] = value
    dirty!(obj)
    return value
end

function set_external_param!(obj::ModelDef, name::Symbol, value::Number;
                             param_dims::Union{Nothing,Array{Symbol}} = nothing)
    set_external_scalar_param!(obj, name, value)
end

function set_external_param!(obj::ModelDef, name::Symbol,
                             value::Union{AbstractArray, AbstractRange, Tuple};
                             param_dims::Union{Nothing,Array{Symbol}} = nothing)
    ti = get_time_index_position(param_dims)
    if ti != nothing
        value = convert(Array{number_type(obj)}, value)
        num_dims = length(param_dims)
        values = get_timestep_array(obj, eltype(value), num_dims, ti, value)
    else
        values = value
    end

    set_external_array_param!(obj, name, values, param_dims)
end

"""
    set_external_array_param!(obj::ModelDef,
                              name::Symbol, value::TimestepVector, dims)

Add a one dimensional time-indexed array parameter indicated by `name` and
`value` to the composite `obj`.  In this case `dims` must be `[:time]`.
"""
function set_external_array_param!(obj::ModelDef,
                                   name::Symbol, value::TimestepVector, dims)
    param = ArrayModelParameter(value, [:time])  # must be :time
    set_external_param!(obj, name, param)
end

"""
    set_external_array_param!(obj::ModelDef,
                              name::Symbol, value::TimestepMatrix, dims)

Add a multi-dimensional time-indexed array parameter `name` with value
`value` to the composite `obj`.  In this case `dims` must be `[:time]`.
"""
function set_external_array_param!(obj::ModelDef,
                                   name::Symbol, value::TimestepArray, dims)
    param = ArrayModelParameter(value, dims === nothing ? Vector{Symbol}() : dims)
    set_external_param!(obj, name, param)
end

"""
    set_external_array_param!(obj::ModelDef,
                              name::Symbol, value::AbstractArray, dims)

Add an array type parameter `name` with value `value` and `dims` dimensions to the composite `obj`.
"""
function set_external_array_param!(obj::ModelDef,
                                   name::Symbol, value::AbstractArray, dims)
    param = ArrayModelParameter(value, dims === nothing ? Vector{Symbol}() : dims)
    set_external_param!(obj, name, param)
end

"""
    set_external_scalar_param!(obj::ModelDef, name::Symbol, value::Any)

Add a scalar type parameter `name` with the value `value` to the composite `obj`.
"""
function set_external_scalar_param!(obj::ModelDef, name::Symbol, value::Any)
    param = ScalarModelParameter(value)
    set_external_param!(obj, name, param)
end

"""
    update_param!(obj::AbstractCompositeComponentDef, name::Symbol, value; update_timesteps = false)

Update the `value` of an external model parameter in composite `obj`, referenced
by `name`. Optional boolean argument `update_timesteps` with default value
`false` indicates whether to update the time keys associated with the parameter
values to match the model's time index.
"""
function update_param!(obj::AbstractCompositeComponentDef, name::Symbol, value; update_timesteps = false)
    _update_param!(obj::AbstractCompositeComponentDef, name, value, update_timesteps; raise_error = true)
end

function _update_param!(obj::AbstractCompositeComponentDef,
                        name::Symbol, value, update_timesteps; raise_error = true)
    param = external_param(obj, name, missing_ok=true)
    if param === nothing
        error("Cannot update parameter; $name not found in composite's external parameters.")
    end

    if param isa ScalarModelParameter
        if update_timesteps && raise_error
            error("Cannot update timesteps; parameter $name is a scalar parameter.")
        end
        _update_scalar_param!(param, name, value)
    else
        _update_array_param!(obj, name, value, update_timesteps, raise_error)
    end

    dirty!(obj)
end

function _update_scalar_param!(param::ScalarModelParameter, name, value)
    if ! (value isa typeof(param.value))
        try
            value = convert(typeof(param.value), value)
        catch e
            error("Cannot update parameter $name; expected type $(typeof(param.value)) but got $(typeof(value)).")
        end
    end
    param.value = value
    nothing
end

function _update_array_param!(obj::AbstractCompositeComponentDef, name, value, update_timesteps, raise_error)
    # Get original parameter
    param = external_param(obj, name)

    # Check type of provided parameter
    if !(typeof(value) <: AbstractArray)
        error("Cannot update array parameter $name with a value of type $(typeof(value)).")

    elseif !(eltype(value) <: eltype(param.values))
        try
            value = convert(Array{eltype(param.values)}, value)
        catch e
            error("Cannot update parameter $name; expected array of type $(eltype(param.values)) but got $(eltype(value)).")
        end
    end

    # Check size of provided parameter
    if update_timesteps && param.values isa TimestepArray
        expected_size = ([length(dim_keys(obj, d)) for d in dim_names(param)]...,)
    else
        expected_size = size(param.values)
    end
    if size(value) != expected_size
        error("Cannot update parameter $name; expected array of size $expected_size but got array of size $(size(value)).")
    end

    if update_timesteps
        if param.values isa TimestepArray
            T = eltype(value)
            N = length(size(value))
            ti = get_time_index_position(param)
            new_timestep_array = get_timestep_array(obj, T, N, ti, value)
            set_external_param!(obj, name, ArrayModelParameter(new_timestep_array, dim_names(param)))

        elseif raise_error
            error("Cannot update timesteps; parameter $name is not a TimestepArray.")
        else
            param.values = value
        end
    else
        if param.values isa TimestepArray
            param.values.data = value
        else
            param.values = value
        end
    end
    dirty!(obj)
    nothing
end

"""
    update_params!(obj::AbstractCompositeComponentDef, parameters::Dict{T, Any};
                   update_timesteps = false) where T

For each (k, v) in the provided `parameters` dictionary, `update_param!`
is called to update the external parameter by name k to value v, with optional
Boolean argument update_timesteps. Each key k must be a symbol or convert to a
symbol matching the name of an external parameter that already exists in the
component definition.
"""
function update_params!(obj::AbstractCompositeComponentDef, parameters::Dict; update_timesteps = false)
    parameters = Dict(Symbol(k) => v for (k, v) in parameters)
    for (param_name, value) in parameters
        _update_param!(obj, param_name, value, update_timesteps; raise_error = false)
    end
    nothing
end

function add_connector_comps!(obj::AbstractCompositeComponentDef)
    conns = internal_param_conns(obj)
    i = 1 # counter to track the number of connector comps added

    for comp_def in compdefs(obj)
        comp_name = nameof(comp_def)
        comp_path = comp_def.comp_path

        # first need to see if we need to add any connector components for this component
        internal_conns  = filter(x -> x.dst_comp_path == comp_path, conns)
        need_conn_comps = filter(x -> x.backup !== nothing, internal_conns)

        # isempty(need_conn_comps) || @info "Need connectors comps: $need_conn_comps"

        for conn in need_conn_comps
            add_backup!(obj, conn.backup)

            num_dims = length(size(external_param(obj, conn.backup).values))

            if ! (num_dims in (1, 2))
                error("Connector components for parameters with > 2 dimensions are not implemented.")
            end

            # Fetch the definition of the appropriate connector commponent
            conn_comp_def = (num_dims == 1 ? Mimi.ConnectorCompVector : Mimi.ConnectorCompMatrix)
            conn_comp_name = connector_comp_name(i) # generate a new name
            i += 1 # increment connector comp counter

            # Add the connector component before the user-defined component that required it
            conn_comp = add_comp!(obj, conn_comp_def, conn_comp_name, before=comp_name)
            conn_path = conn_comp.comp_path

            # add a connection between src_component and the ConnectorComp
            add_internal_param_conn!(obj, InternalParameterConnection(conn.src_comp_path, conn.src_var_name,
                                                                      conn_path, :input1,
                                                                      conn.ignoreunits))

            # add a connection between ConnectorComp and dst_component
            add_internal_param_conn!(obj, InternalParameterConnection(conn_path, :output,
                                                                      conn.dst_comp_path, conn.dst_par_name,
                                                                      conn.ignoreunits))

            # add a connection between ConnectorComp and the external backup data
            add_external_param_conn!(obj, ExternalParameterConnection(conn_path, :input2, conn.backup))

            # TBD: first/last stuff may be deprecated
            src_comp_def = compdef(obj, conn.src_comp_path)
            set_param!(obj, conn_comp_name, :first, first_period(obj, src_comp_def))
            set_param!(obj, conn_comp_name, :last, last_period(obj, src_comp_def))
        end
    end

    # Save the sorted component order for processing
    # obj.sorted_comps = _topological_sort(obj)

    return nothing
end


text:
"""

```
neighborlist(box, cl; parallel=true)
```

Compute the neighbor list of a single set or set pairs of particles. Returns a vector of tuples
with all indices of the particles that are within `box.cutoff`, and the distances.

## Example
```julia-repl
julia> x = [ rand(3) for i in 1:1000 ];

julia> box = Box([1,1,1],0.02)
Box{OrthorhombicCell, 3, Float64, 9}
  unit cell matrix: [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
  cutoff: 0.02
  number of computing cells on each dimension: [52, 52, 52]
  computing cell sizes: [0.02, 0.02, 0.02] (lcell: 1)
  Total number of cells: 140608

julia> cl = CellList(x,box) # single set
CellList{3, Float64}
  999 cells with real particles.
  1126 particles in computing box, including images.

julia> CellListMap.neighborlist(box,cl,parallel=false)
15-element Vector{Tuple{Int64, Int64, Float64}}:
 (187, 511, 0.010346860078531755)
 (203, 708, 0.010777737363239403)
 (296, 579, 0.018124283912224655)
 ⋮
 (584, 4, 0.016935844769524398)
 (725, 749, 0.019971874892397875)
 (773, 119, 0.01835233336121765)
 (927, 8, 0.011234110402648743)

```

To obtain the neighbor list (within the cutoff) between two sets of
particles, initialize the cell lists with the two sets:

```julia-repl
julia> x = [ rand(3) for i in 1:1000 ];

julia> y = [ rand(3) for i in 1:1000 ];

julia> box = Box([1,1,1],0.02);

julia> cl = CellList(x,y,box)

julia> cl = CellList(x,y,box)
CellListMap.CellListPair{Vector{SVector{3, Float64}}, 3, Float64}
   1000 particles in the reference vector.
   997 cells with real particles of target vector.

julia> CellListMap.neighborlist(box,cl)
35-element Vector{Tuple{Int64, Int64, Float64}}:
 (409, 982, 0.01634641594779082)
 (521, 422, 0.00919026348035512)
 (625, 731, 0.012986301890746663)
 ⋮
 (647, 730, 0.01565763971458105)
 (296, 668, 0.016556686306217868)
 (992, 589, 0.018392993428289553)

```

"""
function neighborlist(box::Box, cl; parallel=true)

    # Function adds pair to the list
    function push_pair!(i, j, d2, pairs)
        d = sqrt(d2)
        push!(pairs, (i, j, d))
        return pairs
    end

    # We have to define our own reduce function here (for the parallel version)
    # this reduction can be dum assynchronously on a preallocated array
    function reduce_pairs(pairs, pairs_threaded)
        ranges = cumsum(length.(pairs_threaded))
        npairs = ranges[end]
        pairs = resize!(pairs,npairs)
        @sync for it in 1:length(pairs_threaded)
            range = ranges[it]-length(pairs_threaded[it])+1:ranges[it]
            Threads.@spawn pairs[range] .= pairs_threaded[it]
        end
        return pairs
    end

    # Initialize
    pairs = Tuple{Int,Int,typeof(box.cutoff)}[]

    # Run pairwise computation
    pairs = map_pairwise!(
      (x, y, i, j, d2, pairs) -> push_pair!(i, j, d2, pairs),
      pairs, box, cl,
      reduce=reduce_pairs, parallel=parallel
    )
    return pairs

end

"""

```
neighborlist(x,r;parallel=true)
```

Computes the list of pairs of particles in `x` which are closer to each other than `r`.

## Example
```julia-repl
julia> x = [ rand(3) for i in 1:10_000 ];

julia> CellListMap.neighborlist(x,0.05)
24848-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 1055, 0.022977369806392412)
 (1, 5086, 0.026650609138167428)
 ⋮
 (9989, 3379, 0.0467653507446483)
 (9989, 5935, 0.02432728985151653)

```

"""
function neighborlist(x,r;parallel=true)
    box = Box(limits(x),r)
    cl = CellList(x,box,parallel=parallel)
    return neighborlist(box,cl,parallel=parallel)
end

"""

```
neighborlist(x,y,r;parallel=true,autoswap=true)
```

Computes the list of pairs of particles of `x` which are closer than `r` to
the particles of `y`. The `autoswap` option will swap `x` and `y` to try to optimize
the cost of the construction of the cell list.

## Example
```julia-repl
julia> x = [ rand(3) for i in 1:10_000 ];

julia> y = [ rand(3) for i in 1:1_000 ];

julia> CellListMap.neighborlist(x,y,0.05)
5006-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 269, 0.04770884036497686)
 (25, 892, 0.03850515231540869)
 ⋮
 (9952, 749, 0.048875643578313456)
 (9984, 620, 0.04101242499363183)

```

"""
function neighborlist(x,y,r;parallel=true,autoswap=true)
    box = Box(limits(x,y),r)
    cl = CellList(x,y,box,parallel=parallel,autoswap=autoswap)
    return neighborlist(box,cl,parallel=parallel)
end


text:
import Base: convert

export restorecontext

"""
    restorecontext(term::DeBruijn.Term)

Construct a `NamingContext` with made-up names for the free variables in `term` (respecting their
actual indices).
"""
function restorecontext(t::DeBruijn.Term; freenamehint = :free, boundnamehint = :x)
    @assert freenamehint != boundnamehint
    m = reduce(max, freevars(t), init = 0)
    NamingContext((Symbol(freenamehint, i) for i = 1:m), namehint = boundnamehint)
end

# Named -> DeBruijn
convert(::Type{DeBruijn.Term}, t::Named.Term) =
    convert(DeBruijn.Term, t, NamingContext(freevars(t)))
convert(::Type{DeBruijn.Term}, v::Named.Var, Γ::NamingContext) =
    let index = findfirst(isequal(v.name), Γ)
        index !== nothing ? DeBruijn.Var(index) : throw(KeyError(v.name))
    end
convert(::Type{DeBruijn.Term}, t::Named.App, Γ::NamingContext) =
    let (Γₗ, Γᵣ) = split(Γ)
        DeBruijn.App(convert(DeBruijn.Term, t.car, Γₗ), convert(DeBruijn.Term, t.cdr, Γᵣ))
    end
convert(::Type{DeBruijn.Term}, t::Named.Abs, Γ::NamingContext) =
        DeBruijn.Abs(convert(DeBruijn.Term, t.body, pushfirst(Γ, t.boundname)))

# Named -> LocallyNameless
convert(::Type{LocallyNameless.Term}, t::Named.Term) =
    convert(LocallyNameless.Term, t, NamingContext(freevars(t)))
convert(::Type{LocallyNameless.Term}, v::Named.Var, Γ::NamingContext, level::Int = 0) =
    let index = findfirst(isequal(v.name), Γ)
        index > level ? LocallyNameless.FVar(Γ[index]) : LocallyNameless.BVar(index)
    end
convert(::Type{LocallyNameless.Term}, t::Named.App, Γ::NamingContext, level::Int = 0) =
    let (Γₗ, Γᵣ) = split(Γ)
        LocallyNameless.App(convert(LocallyNameless.Term, t.car, Γₗ, level),
                            convert(LocallyNameless.Term, t.cdr, Γᵣ, level))
    end
convert(::Type{LocallyNameless.Term}, t::Named.Abs, Γ::NamingContext, level::Int = 0) =
    LocallyNameless.Abs(convert(LocallyNameless.Term, t.body, pushfirst(Γ, t.boundname), level + 1))

# DeBruijn -> Named
convert(::Type{Named.Term}, t::DeBruijn.Term) = convert(Named.Term, t, restorecontext(t))
convert(::Type{Named.Term}, t::DeBruijn.Var, Γ::NamingContext) = Named.Var(Γ[t.index])
convert(::Type{Named.Term}, t::DeBruijn.App, Γ::NamingContext) =
    let (Γₗ, Γᵣ) = split(Γ)
        Named.App(convert(Named.Term, t.car, Γₗ), convert(Named.Term, t.cdr, Γᵣ))
    end
convert(::Type{Named.Term}, t::DeBruijn.Abs, Γ::NamingContext) =
    let (boundname, Γ′) = freshname(Γ)
        Named.Abs(boundname, convert(Named.Term, t.body, Γ′))
    end

# DeBruijn -> LocallyNameless
convert(::Type{LocallyNameless.Term}, t::DeBruijn.Term) =
    convert(LocallyNameless.Term, t, restorecontext(t))
convert(::Type{LocallyNameless.Term}, t::DeBruijn.Var, Γ::NamingContext, level::Int = 0) =
    t.index > level ? LocallyNameless.FVar(Γ[t.index - level]) : LocallyNameless.BVar(t.index)
convert(::Type{LocallyNameless.Term}, t::DeBruijn.App, Γ::NamingContext, level::Int = 0) =
    LocallyNameless.App(convert(LocallyNameless.Term, t.car, Γ, level),
                        convert(LocallyNameless.Term, t.cdr, Γ, level))
convert(::Type{LocallyNameless.Term}, t::DeBruijn.Abs, Γ::NamingContext, level::Int = 0) =
    LocallyNameless.Abs(convert(LocallyNameless.Term, t.body, Γ, level + 1))

# LocallyNameless -> Named
convert(::Type{Named.Term}, t::LocallyNameless.Term) =
    convert(Named.Term, t, NamingContext(freevars(t)))
convert(::Type{Named.Term}, t::LocallyNameless.FVar, Γ::NamingContext) = Named.Var(t.name)
convert(::Type{Named.Term}, t::LocallyNameless.BVar, Γ::NamingContext) = Named.Var(Γ[t.index])
convert(::Type{Named.Term}, t::LocallyNameless.App, Γ::NamingContext) =
    let (Γₗ, Γᵣ) = split(Γ)
        Named.App(convert(Named.Term, t.car, Γₗ), convert(Named.Term, t.cdr, Γᵣ))
    end
convert(::Type{Named.Term}, t::LocallyNameless.Abs, Γ::NamingContext) =
    let (boundname, Γ′) = freshname(Γ)
        Named.Abs(boundname, convert(Named.Term, t.body, Γ′))
    end

# LocallyNameless -> DeBruijn
convert(::Type{DeBruijn.Term}, t::LocallyNameless.Term) =
    convert(DeBruijn.Term, t, NamingContext(freevars(t)))
convert(::Type{DeBruijn.Term}, v::LocallyNameless.BVar, Γ::NamingContext) = DeBruijn.Var(v.index)
convert(::Type{DeBruijn.Term}, v::LocallyNameless.FVar, Γ::NamingContext) =
    let index = findfirst(isequal(v.name), Γ)
        index !== nothing ? DeBruijn.Var(index) : throw(KeyError(v.name))
    end
convert(::Type{DeBruijn.Term}, t::LocallyNameless.App, Γ::NamingContext) =
    let (Γₗ, Γᵣ) = split(Γ)
        DeBruijn.App(convert(DeBruijn.Term, t.car, Γₗ), convert(DeBruijn.Term, t.cdr, Γᵣ))
    end
convert(::Type{DeBruijn.Term}, t::LocallyNameless.Abs, Γ::NamingContext) =
    let (boundname, Γ′) = freshname(Γ)
        DeBruijn.Abs(convert(DeBruijn.Term, t.body, Γ′))
    end