text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/optimize/optimize_convex.jl
# Contains the solve_local_nlp! routine which computes the optimal value
# of a convex function. This is used to compute the upper bound in the
# branch and bound routine. A number of utility function required for
# solve_local_nlp! are also included.
#############################################################################

function _update_branch_variables!(nlp_opt, m)
    n = m._current_node
    for i = 1:m._working_problem._variable_num
        vinfo = @inbounds m._working_problem._variable_info[i]
        if m._branch_variables[i]
            if vinfo.is_integer
            else
                indx = @inbounds m._sol_to_branch_map[i]
                lvb  = @inbounds n.lower_variable_bound[indx]
                uvb  = @inbounds n.upper_variable_bound[indx]
                if vinfo.is_fixed
                    MOI.add_constraint(nlp_opt, SV(m.upper_variables[i]), ET(lvb))
                elseif vinfo.has_lower_bound
                    if vinfo.has_upper_bound
                        MOI.add_constraint(nlp_opt, SV(m.upper_variables[i]), LT(uvb))
                        MOI.add_constraint(nlp_opt, SV(m.upper_variables[i]), GT(lvb))
                    else
                        MOI.add_constraint(nlp_opt, SV(m.upper_variables[i]), GT(lvb))
                    end
                elseif vinfo.has_upper_bound
                    MOI.add_constraint(nlp_opt, SV(m.upper_variables[i]), LT(uvb))
                end
            end
        end
    end
    return nothing
end

# TODO: Select between infeasible starting point reduction strategies
function _set_starting_point!(opt, m)
    for i = 1:m._working_problem._variable_num
        if m._branch_variables[i]
            vp = mid(m, m._sol_to_branch_map[i])
        else
            vp = mid(m._working_problem.variable_info[i])
            if isinf(vp)
                if isinf(_lower_bound(m._working_problem.variable_info[i]))
                    vp = _lower_bound(m._working_problem.variable_info[i])
                else
                    vp = _upper_bound(m._working_problem.variable_info[i])
                end
            end
        end
        MOI.set(opt, MOI.VariablePrimalStart(m._upper_variables[i]), vp)
    end
    return nothing
end

"""
$(SIGNATURES)

Takes an `MOI.TerminationStatusCode` and a `MOI.ResultStatusCode` and returns `true`
if this corresponds to a solution that is proven to be feasible.
Returns `false` otherwise.
"""
function is_feasible_solution(t::MOI.TerminationStatusCode, r::MOI.ResultStatusCode)

    termination_flag = false
    result_flag = false

    (t == MOI.OPTIMAL) && (termination_flag = true)
    (t == MOI.LOCALLY_SOLVED) && (termination_flag = true)

    # This is default solver specific... the acceptable constraint tolerances
    # are set to the same values as the basic tolerance. As a result, an
    # acceptably solved solution is feasible but non necessarily optimal
    # so it should be treated as a feasible point
    if (t == MOI.ALMOST_LOCALLY_SOLVED) && (r == MOI.NEARLY_FEASIBLE_POINT)
        termination_flag = true
        result_flag = true
    end

    (r == MOI.FEASIBLE_POINT) && (result_flag = true)

    return (termination_flag && result_flag)
end

function _unpack_local_nlp_solve!(m::Optimizer, opt::T, idx_map; adjust_bnd::Bool = true) where T
    # Process output info and save to CurrentUpperInfo object
    m._upper_termination_status = MOI.get(nlp_optimizer, MOI.TerminationStatus())
    m._upper_result_status = MOI.get(nlp_optimizer, MOI.PrimalStatus())
    if is_feasible_solution(m._upper_termination_status, m._upper_result_status)
        m._upper_feasibility = true
        value = MOI.get(nlp_optimizer, MOI.ObjectiveValue())
        m._best_upper_value = min(value, m._best_upper_value)
        m._upper_solution .= MOI.get(nlp_optimizer, MOI.VariablePrimal(), upper_variables)
    else
        m._upper_feasibility = false
        m._upper_objective_value = Inf
    end
    return nothing
end

"""

Constructs and solves the problem locally on on node `y` updated the upper
solution informaton in the optimizer.
"""
function solve_local_nlp!(m::Optimizer, i)

    opt = m.nlp_optimizer
    #set_config!(m, opt)
    bridged_opt = _bridge_optimizer(Val{DIFF_CVX}(), opt)
    idx_map = MOIU.default_copy_to(bridged_opt, m._input_problem, false)
    (m.verbosity < 5) && MOI.set(bridged_opt, MOI.Silent(), true)

    _update_branch_variables!(bridged_opt, m, i)
    _set_starting_point!(bridged_opt, m, i)
    MOI.optimize!(bridged_opt)
    _unpack_local_nlp_solve!(m, bridged_opt, idx_map)

    return nothing
end

"""
$(SIGNATURES)

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

"""
$(SIGNATURES)

Default upper bounding problem which simply calls `solve_local_nlp!` to solve
the nlp locally.
"""
function upper_problem!(t::ExtensionType, m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}

    if !default_nlp_heurestic(m)
        m._upper_feasibility = false
        m._upper_objective_value = Inf
    elseif _multistart_heurestic(m)
        _multistart_local_nlp!(m)
    else
        _solve_local_nlp!(m, 1)
    end
    return nothing
end


text:
###################
# Pretty Printing #
###################

function prettytime(t)
    if t < 1e3
        value, units = t, "ns"
    elseif t < 1e6
        value, units = t / 1e3, "μs"
    elseif t < 1e9
        value, units = t / 1e6, "ms"
    else
        value, units = t / 1e9, "s"
    end

    if round(value) >= 100
        str = string(@sprintf("%.0f", value), units)
    elseif round(value * 10) >= 100
        str = string(@sprintf("%.1f", value), units)
    else
        str = string(@sprintf("%.2f", value), units)
    end
    return lpad(str, 6, " ")
end

function prettymemory(b)
    if b < 1000
        value, units = b, "B"
    elseif b < 1000^2
        value, units = b / 1024, "KiB"
    elseif b < 1000^3
        value, units = b / 1024^2, "MiB"
    else
        value, units = b / 1024^3, "GiB"
    end

    if round(value) >= 100
        str = string(@sprintf("%.0f", value), units)
    elseif round(value * 10) >= 100
        str = string(@sprintf("%.1f", value), units)
    elseif value >= 0
        str = string(@sprintf("%.2f", value), units)
    else
        str = "-"
    end
    return lpad(str, 7, " ")
end

function prettypercent(nominator, denominator)
    value = nominator / denominator * 100

    if denominator == 0 && nominator == 0
        str = " - %"
    elseif denominator == 0
        str = "inf %"
    elseif round(value) >= 100
        str = string(@sprintf("%.0f", value), "%")
    elseif round(value * 10) >= 100
        str = string(@sprintf("%.1f", value), "%")
    else
        str = string(@sprintf("%.2f", value), "%")
    end
    return str
end

function prettycount(t::Integer)
    if t < 1000
        return string(t)
    elseif t < 1000^2
        value, units = t / 1000, "k"
    elseif t < 1000^3
        value, units = t / 1e6, "M"
    else
        value, units = t / 1e9, "B"
    end

    if round(value) >= 100
        str = string(@sprintf("%.0f", value), units)
    elseif round(value * 10) >= 100
        str = string(@sprintf("%.1f", value), units)
    else
        str = string(@sprintf("%.2f", value), units)
    end
    return str
end

function rpad(
    s::Union{AbstractChar,AbstractString},
    n::Integer,
    p::Union{AbstractChar,AbstractString}=' ',
) :: String
    n = Int(n)::Int
    m = signed(n) - Int(textwidth(s))::Int
    m ≤ 0 && return string(s)
    l = textwidth(p)
    q, r = divrem(m, l)
    r == 0 ? string(s, p^q) : string(s, p^q, first(p, r))
end

ext: jl


text:
module Distances

using LinearAlgebra
using Statistics
import StatsAPI: pairwise, pairwise!

export
    # generic types/functions
    PreMetric,
    SemiMetric,
    Metric,

    # generic functions
    result_type,
    colwise,
    pairwise,
    colwise!,
    pairwise!,
    evaluate,

    # distance classes
    Euclidean,
    SqEuclidean,
    PeriodicEuclidean,
    Cityblock,
    TotalVariation,
    Chebyshev,
    Minkowski,
    Jaccard,
    BrayCurtis,
    RogersTanimoto,

    Hamming,
    CosineDist,
    CorrDist,
    ChiSqDist,
    KLDivergence,
    GenKLDivergence,
    JSDivergence,
    RenyiDivergence,
    SpanNormDist,

    WeightedEuclidean,
    WeightedSqEuclidean,
    WeightedCityblock,
    WeightedMinkowski,
    WeightedHamming,
    SqMahalanobis,
    Mahalanobis,
    BhattacharyyaDist,
    HellingerDist,

    Haversine,
    SphericalAngle,

    MeanAbsDeviation,
    MeanSqDeviation,
    RMSDeviation,
    NormRMSDeviation,
    Bregman,

    # convenient functions
    euclidean,
    sqeuclidean,
    peuclidean,
    cityblock,
    totalvariation,
    jaccard,
    braycurtis,
    rogerstanimoto,
    chebyshev,
    minkowski,

    hamming,
    cosine_dist,
    corr_dist,
    chisq_dist,
    kl_divergence,
    gkl_divergence,
    js_divergence,
    renyi_divergence,
    spannorm_dist,

    weuclidean,
    wsqeuclidean,
    wcityblock,
    wminkowski,
    whamming,
    sqmahalanobis,
    mahalanobis,
    bhattacharyya,
    hellinger,
    bregman,

    haversine,
    spherical_angle,

    meanad,
    msd,
    rmsd,
    nrmsd

if VERSION < v"1.2-"
    import Base: has_offset_axes
    require_one_based_indexing(A...) =
        !has_offset_axes(A...) ||
            throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
else
    import Base: require_one_based_indexing
end

include("common.jl")
include("generic.jl")
include("metrics.jl")
include("haversine.jl")
include("mahalanobis.jl")
include("bhattacharyya.jl")
include("bregman.jl")

end # module end


text:
"Compilation to Julia (and other targets, if you want)"
module Compile

using ..AExpressions, ..CompileUtils
import MacroTools: striplines

export compiletojulia, runprogram

"compile `aexpr` into Expr"
function compiletojulia(aexpr::AExpr)::Expr

  # dictionary containing types/definitions of global variables, for use in constructing init func.,
  # next func., etcetera; the three categories of global variable are external, initnext, and lifted
  historydata = Dict([("external" => []), # :typedecl aexprs for all external variables
               ("initnext" => []), # :assign aexprs for all initnext variables
               ("lifted" => []), # :assign aexprs for all lifted variables
               ("types" => Dict())]) # map of global variable names (symbols) to types

  if (aexpr.head == :program)
    # handle AExpression lines
    lines = filter(x -> x !== :(), map(arg -> compile(arg, historydata, aexpr), aexpr.args))

    # construct STATE struct and initialize state::STATE
    stateStruct = compilestatestruct(historydata)
    initStateStruct = compileinitstate(historydata)

    # handle init, next, prev, and built-in functions
    initnextFunctions = compileinitnext(historydata)
    prevFunctions = compileprevfuncs(historydata)
    builtinFunctions = compilebuiltin()

    # remove empty lines
    lines = filter(x -> x != :(),
            vcat(builtinFunctions, lines, stateStruct, initStateStruct, prevFunctions, initnextFunctions))

    # construct module
    expr = quote
      module CompiledProgram
        export init, next
        import Base.min
        using Distributions
        using MLStyle: @match
        $(lines...)
      end
    end
    expr.head = :toplevel
    striplines(expr)
  else
    throw(AutumnCompileError("AExpr Head != :program"))
  end
end

"Run `prog` for finite number of time steps"
function runprogram(prog::Expr, n::Int)
  mod = eval(prog)
  state = mod.init(mod.Click(5, 5))

  for i in 1:n
    externals = [nothing, mod.Click(rand([1:10;]), rand([1:10;]))]
    state = mod.next(mod.next(state, externals[rand(Categorical([0.7, 0.3]))]))
  end
  state
end

end


text:
export ref_count

import Base: show

"""
    ref_count()

Make a ConnectableObservable behave like a ordinary observable and automates the way you can connect to it.
Internally it counts the subscriptions to the observable and subscribes (only once) to the source if the number of subscriptions is larger than 0. If the number of subscriptions is smaller than 1, it unsubscribes from the source.
This way you can make sure that everything before the published refCount has only a single subscription independently of the number of subscribers to the target observable.

Note that using the [`share`](@ref) operator is exactly the same as using the [`publish`](@ref) operator (making the observable hot) and the [`ref_count()`](@ref) operator in a sequence.

# Example

```jldoctest
using Rocket

subject = Subject(Int, scheduler = AsapScheduler())
source  = from(1:5) |> multicast(subject) |> ref_count()

actor1 = logger("1")
actor2 = logger("2")

subscription1 = subscribe!(source, actor1)
subscription2 = subscribe!(source, actor2)

unsubscribe!(subscription1)
unsubscribe!(subscription2)
;

# output
[1] Data: 1
[1] Data: 2
[1] Data: 3
[1] Data: 4
[1] Data: 5
[1] Completed
[2] Completed
```

See also: [`AbstractOperator`](@ref), [`publish`](@ref), [`multicast`](@ref), [`share`](@ref)
"""
ref_count() = RefCountOperator()

struct RefCountOperator <: InferableOperator end

function on_call!(::Type{L}, ::Type{L}, operator::RefCountOperator, source) where L
    return proxy(L, source, RefCountProxy())
end

operator_right(operator::RefCountOperator, ::Type{L}) where L = L

struct RefCountProxy <: SourceProxy end

source_proxy!(::Type{L}, proxy::RefCountProxy, source::S) where { L, S } = RefCountSource{L, S}(source)

@subscribable mutable struct RefCountSource{L, S} <: Subscribable{L}
    csource       :: S
    refcount      :: Int
    csubscription :: Teardown

    RefCountSource{L, S}(csource::S) where { L, S } = new(csource, 0, voidTeardown)
end

function on_subscribe!(refcount_source::RefCountSource, actor)
    subscription = subscribe!(refcount_source.csource, actor)
    refcount_source.refcount += 1
    if refcount_source.refcount === 1
        refcount_source.csubscription = connect(refcount_source.csource)
    end
    return RefCountSourceSubscription(refcount_source, subscription)
end

mutable struct RefCountSourceSubscription{S, T} <: Teardown
    refcount_source :: Union{Nothing, S}
    subscription    :: T
end

as_teardown(::Type{<:RefCountSourceSubscription}) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::RefCountSourceSubscription)
    unsubscribe!(subscription.subscription)

    refcount_source = subscription.refcount_source
    if refcount_source !== nothing
        refcount_source.refcount -= 1
        if refcount_source.refcount === 0
            unsubscribe!(refcount_source.csubscription)
            refcount_source.csubscription = voidTeardown
        end
        subscription.refcount_source = nothing
    end
end

Base.show(io::IO, ::RefCountOperator)           = print(io, "RefCountOperator()")
Base.show(io::IO, ::RefCountProxy)              = print(io, "RefCountProxy()")
Base.show(io::IO, ::RefCountSource{L})  where L = print(io, "RefCountSource($L)")
Base.show(io::IO, ::RefCountSourceSubscription) = print(io, "RefCountSubscription()")


text:
module IndexSetsTests

using PartitionedArrays
using Test

part = 2
noids = 20
firstgid = 11
ids = IndexRange(part,noids,firstgid)
@test ids.part == part
@test ids.lid_to_gid == 11:30
push!(ids.lid_to_gid,48)
push!(ids.lid_to_gid,49)
@test ids.lid_to_gid == vcat(collect(11:30),[48,49])
@test ids.lid_to_part == fill(part,20)
push!(ids.lid_to_part,3)
push!(ids.lid_to_part,3)
@test ids.lid_to_part == vcat(fill(part,20),[3,3])
@test ids.oid_to_lid == 1:20
@test ids.hid_to_lid == Int32[]
push!(ids.hid_to_lid,21)
push!(ids.hid_to_lid,22)
@test ids.hid_to_lid == Int32[21,22]
@test ids.lid_to_ohid == 1:20
push!(ids.lid_to_ohid,-1)
push!(ids.lid_to_ohid,-2)
@test ids.lid_to_ohid == vcat(collect(1:20),[-1,-2])
@test collect(keys(ids.gid_to_lid)) == 11:30
@test collect(values(ids.gid_to_lid)) == 1:20
ids.gid_to_lid[48] = 21
ids.gid_to_lid[49] = 22
@test haskey(ids.gid_to_lid,49)
@test haskey(ids.gid_to_lid,48)
@test haskey(ids.gid_to_lid,12)
@test haskey(ids.gid_to_lid,30)
@test haskey(ids.gid_to_lid,11)
@test !haskey(ids.gid_to_lid,10)
@test !haskey(ids.gid_to_lid,31)
@test length(collect(ids.gid_to_lid)) == 22

hid_to_gid = [48,49]
hid_to_part = Int32[3,3]
ids = IndexRange(part,noids,firstgid,hid_to_gid,hid_to_part)
@test ids.part == part
@test ids.lid_to_gid == vcat(collect(11:30),[48,49])
@test ids.lid_to_part == vcat(fill(part,20),[3,3])
@test ids.hid_to_lid == Int32[21,22]
@test ids.lid_to_ohid == vcat(collect(1:20),[-1,-2])
@test haskey(ids.gid_to_lid,49)
@test haskey(ids.gid_to_lid,48)
@test haskey(ids.gid_to_lid,12)
@test haskey(ids.gid_to_lid,30)
@test haskey(ids.gid_to_lid,11)
@test !haskey(ids.gid_to_lid,10)
@test !haskey(ids.gid_to_lid,31)
@test length(collect(ids.gid_to_lid)) == 22

lid_to_gid = collect(ids.lid_to_gid)
lid_to_part = collect(ids.lid_to_part)
ids = ExtendedIndexRange(part,lid_to_gid,lid_to_part,firstgid)
@test ids.part == part
@test ids.lid_to_gid == vcat(collect(11:30),[48,49])
@test ids.lid_to_part == vcat(fill(part,20),[3,3])
@test ids.hid_to_lid == Int32[21,22]
@test ids.lid_to_ohid == vcat(collect(1:20),[-1,-2])
@test haskey(ids.gid_to_lid,49)
@test haskey(ids.gid_to_lid,48)
@test haskey(ids.gid_to_lid,12)
@test haskey(ids.gid_to_lid,30)
@test haskey(ids.gid_to_lid,11)
@test !haskey(ids.gid_to_lid,10)
@test !haskey(ids.gid_to_lid,31)
@test length(collect(ids.gid_to_lid)) == 22
push!(ids.lid_to_gid,62)
push!(ids.lid_to_gid,63)
push!(ids.lid_to_part,4)
push!(ids.lid_to_part,4)
push!(ids.hid_to_lid,23)
push!(ids.hid_to_lid,24)
push!(ids.lid_to_ohid,-3)
push!(ids.lid_to_ohid,-4)
ids.gid_to_lid[62] = 23
ids.gid_to_lid[63] = 24
@test ids.lid_to_gid == vcat(collect(11:30),[48,49,62,63])
@test ids.lid_to_part == vcat(fill(part,20),[3,3,4,4])
@test ids.hid_to_lid == Int32[21,22,23,24]
@test ids.lid_to_ohid == vcat(collect(1:20),[-1,-2,-3,-4])
@test haskey(ids.gid_to_lid,62)
@test haskey(ids.gid_to_lid,63)
@test !haskey(ids.gid_to_lid,64)
@test !haskey(ids.gid_to_lid,61)
@test length(collect(ids.gid_to_lid)) == 24

end # module


text:
""" Categories of graphs and other categorical and algebraic aspects of graphs.
"""
module GraphCategories

using DataStructures

using ..FinSets, ...ACSetInterface, ..Limits
using ...Graphs.BasicGraphs
import ...Graphs.GraphAlgorithms: connected_component_projection,
  connected_component_projection_bfs

function connected_component_projection(g::ACSet)::FinFunction
  proj(coequalizer(FinFunction(src(g), nv(g)),
                   FinFunction(tgt(g), nv(g))))
end

# This algorithm is linear in the number of vertices of g, so it should be
# significantly faster than the previous one in some cases.
function connected_component_projection_bfs(g::ACSet)
  label = zeros(Int, nv(g))

  q = Queue{Int}()
  for v in 1:nv(g)
    label[v] != 0 && continue
    label[v] = v
    empty!(q)
    enqueue!(q, v)
    while !isempty(q)
      src = dequeue!(q)
      for vertex in neighbors(g, src)
        if label[vertex] == 0
          enqueue!(q,vertex)
          label[vertex] = v
        end
      end
    end
  end

  normalized = searchsortedfirst.(Ref(unique!(sort(label))), label)
  FinFunction(normalized)
end

end

max_forks_repo_name: epatters/CompCat


text:
using Optim, Zygote, ParameterHandling

# L-BFGS parameters chosen because they seem to work well empirically.
# You could also try with the defaults.
default_optimizer = LBFGS(;
    alphaguess=Optim.LineSearches.InitialStatic(; scaled=true),
    linesearch=Optim.LineSearches.BackTracking(),
)

function optimize_loss(loss, θ_init; optimizer=default_optimizer, maxiter=1_000)
    options = Optim.Options(; iterations=maxiter, show_trace=true)

    θ_flat_init, unflatten = ParameterHandling.value_flatten(θ_init)
    loss_packed = loss ∘ unflatten

    # https://julianlsolvers.github.io/Optim.jl/stable/#user/tipsandtricks/#avoid-repeating-computations
    function fg!(F, G, x)
        if F != nothing && G != nothing
            val, grad = Zygote.withgradient(loss_packed, x)
            G .= only(grad)
            return val
        elseif G != nothing
            grad = Zygote.gradient(loss_packed, x)
            G .= only(grad)
            return nothing
        elseif F != nothing
            return loss_packed(x)
          end
    end

    result = optimize(
        Optim.only_fg!(fg!),
        θ_flat_init,
        optimizer,
        options;
        inplace=false,
    )

    return unflatten(result.minimizer), result
end


text:
using Kosmos
using Test

@testset "Kosmos.jl" begin

    include("clock.jl")
    include("sysgeneration.jl")

end


text:
# Index manipulations
#---------------------
"""
    permuteind(tsrc::AbstractTensorMap{S}, p1::NTuple{N₁,Int}, p2::NTuple{N₂,Int} = ()) -> tdst::TensorMap{S,N₁,N₂}

Permute the indices of `tsrc::AbstractTensorMap{S}` such that a new tensor
`tdst::TensorMap{S,N₁,N₂}` is obtained, with indices in `p1` playing the role of
the codomain or range of the map, and indices in `p2` indicating the domain.

To permute into an existing `tdst`, use `permuteind!(tdst, tsrc, p1, p2)`.
"""
function permuteind(t::TensorMap{S},
                    p1::IndexTuple{N₁},  p2::IndexTuple{N₂}=();
                    copy::Bool = false) where {S,N₁,N₂}
    cod = ProductSpace{S,N₁}(map(n->space(t, n), p1))
    dom = ProductSpace{S,N₂}(map(n->dual(space(t, n)), p2))

    if !copy
        # share data if possible
        if p1 === codomainind(t) && p2 === domainind(t)
            return t
        elseif has_shared_permuteind(t, p1, p2)
            return TensorMap(reshape(t.data, dim(cod), dim(dom)), cod, dom)
        end
    end
    # general case
    @inbounds begin
        return permuteind!(similar(t, cod←dom), t, p1, p2)
    end
end

function permuteind(t::AdjointTensorMap{S}, p1::IndexTuple{N₁}, p2::IndexTuple{N₂}=();
                    copy::Bool = false) where {S,N₁,N₂}
    p1′ = map(n->adjointtensorindex(t, n), p2)
    p2′ = map(n->adjointtensorindex(t, n), p1)
    adjoint(permuteind(adjoint(t), p1′, p2′; copy = copy))
end

function has_shared_permuteind(t::TensorMap, p1, p2)
    if p1 === codomainind(t) && p2 === domainind(t)
        return true
    elseif sectortype(t) === Trivial
        stridet = i->stride(t[], i)
        sizet = i->size(t[], i)
        canfuse1, d1, s1 = TO._canfuse(sizet.(p1), stridet.(p1))
        canfuse2, d2, s2 = TO._canfuse(sizet.(p2), stridet.(p2))
        return canfuse1 && canfuse2 && s1 == 1 && (d2 == 1 || s2 == d1)
    else
        return false
    end
end

function has_shared_permuteind(t::AdjointTensorMap, p1, p2)
    p1′ = adjointtensorindices(t, p2)
    p2′ = adjointtensorindices(t, p1)
    return has_shared_permuteind(t', p1′, p2′)
end

@propagate_inbounds permuteind!(tdst::AbstractTensorMap{S,N₁,N₂},
                                tsrc::AbstractTensorMap{S},
                                p1::IndexTuple{N₁},
                                p2::IndexTuple{N₂}=()) where {S,N₁,N₂} =
                                add!(true, tsrc, false, tdst, p1, p2)

# Index manipulations that increase or decrease the number of indices

function splitind end#

function fuseind end

# # TODO: reconsider whether we need repartitionind!, or if we just want permuteind!
# function repartitionind!(tdst::TensorMap{S,N₁,N₂}, tsrc::TensorMap{S,N₁′,N₂′}) where {S,N₁,N₂,N₁′,N₂′}
#     space1 = codomain(tdst) ⊗ dual(domain(tdst))
#     space2 = codomain(tsrc) ⊗ dual(domain(tsrc))
#     space1 == space2 || throw(SpaceMismatch())
#     p = (ntuple(n->n, StaticLength(N₁′))..., ntuple(n->N₁′+N₂′+1-n, StaticLength(N₂′)))
#     p1 = TupleTools.getindices(p, ntuple(n->n, StaticLength(N₁)))
#     p2 = reverse(TupleTools.getindices(p, ntuple(n->N₁+n, StaticLength(N₂))))
#     pdata = (p1..., p2...)
#
#     if sectortype(S) == Trivial
#         TensorOperations.add!(1, tsrc[], Val{:N}, 0, tdst[], pdata)
#     else
#         fill!(tdst, 0)
#         for (f1,f2) in fusiontrees(t)
#             for ((f1′,f2′), coeff) in repartition(f1, f2, StaticLength(N₁))
#                 TensorOperations.add!(coeff, tsrc[f1,f2], Val{:N}, 1, tdst[f1′,f2′], pdata)
#             end
#         end
#     end
#     return tdst
# end


text:
module EmpiricalModeDecomposition

# package code goes here
using DSP
using Images
using Interpolations
using Dierckx

include("util.jl")
include("IntrinsicModeFunctions.jl")
include("InstantaneousFrequency.jl")

export findExtrema
export IMF
export IF

end # module

max_stars_repo_name: alexmorley/EmpiricalModeDecomposition

max_issues_repo_name: alexmorley/EmpiricalModeDecomposition

max_forks_repo_name: alexmorley/EmpiricalModeDecomposition


text:
import SExpressions.Keywords: Keyword

@testset "Parser" begin

sx = SExpressions.parse

@testset "numbers" begin
    @test sx("-1") == -1
    @test sx("-1") isa BigInt
    @test sx("+1000") == 1000
    @test sx("+1000") isa BigInt
    @test sx("-15/5") == -3
    @test sx("-15/5") isa Rational{BigInt}
    @test sx("+13/17") == 13//17
    @test sx("+13/17") isa Rational{BigInt}
    @test sx("-1/2") == -1//2
    @test sx("-1/2") isa Rational{BigInt}
    @test sx("-1/2+3i") == -1//2+3im
    @test sx("-1/2+3i") isa Complex{Rational{BigInt}}
    @test sx("3+4/3i") == 3+4im//3
    @test sx("3+4/3i") isa Complex{Rational{BigInt}}
    @test sx("3-2i") == 3-2im
    @test sx("3-2i") isa Complex{BigInt}
    @test_broken sx("0.5") === 0.5
    @test_broken sx("#e0.5") == 1//2
    @test_broken sx("#e0.5") isa Rational{BigInt}
    @test_broken sx("#x03BB") == 955
    @test_broken sx("#x03BB") isa BigInt
end

@testset "booleans" begin
    @test sx("#t")
    @test sx("#true")
    @test sx("#T")
    @test !sx("#f")
    @test !sx("#false")
    @test !sx("#F")
end

@testset "arithmetic" begin
    @test sx("(+ 1 1)") == List(:(+), 1, 1)
    @test sx("(+ 1/2 1/3)") == List(:(+), 1//2, 1//3)
    @test sx("""
    (define (sqr x) (^ x 2))
    """) == SExpression((
            :define,
            (:sqr, :x),
            (:(^), :x, 2)))
end

@testset "invalid syntax" begin
    @test_throws ErrorException SExpressions.parseall("(+ 1 1")
    @test_throws ErrorException SExpressions.parseall("(+ 1 1]")
    @test_throws ErrorException SExpressions.parseall("(+ 1 1))")
    @test_throws ErrorException SExpressions.parseall("(+ 1 1;)")
end

@testset "comment" begin
    @test sx("1;") == 1
    @test sx(";
             x") == :x
    @test sx(";
    ;
    ;
    4;2
    ") == 4
end

@testset "strings" begin
    @test sx("\"x\"") == "x"
    @test sx("\"x\\n\"") == "x\n"
    @test sx("""\"
    x\\
    y\\\\
    z\"""") == "\nxy\\\nz"
end

@testset "symbols" begin
    @test sx("foo") == :foo
    @test_broken sx("|x|") == :x
    @test_broken sx("| |") == Symbol(" ")
end

@testset "lists" begin
    @test sx("(1 2 3)") == List(1, 2, 3)
    @test sx("(1 . 2)") == Cons(1, 2)
    @test sx("(1 . (2))") == List(1, 2)
    @test sx("(1 . + . 2)") == List(:+, 1, 2)
    @test sx("(1 2 . 3)") == Cons(1, Cons(2, 3))
    @test_throws ErrorException sx("(.)")
    @test_throws ErrorException sx("(. x)")
    @test_throws ErrorException sx("(1 . 2 3)")
    @test_throws ErrorException sx("(1 . 2 . 3 . 4)")
end

@testset "reader macros" begin
    @test sx("'x") == List(:quote, :x)
    @test sx("`(+ 1 2)") == List(:quasiquote, List(:(+), 1, 2))
    @test sx("`(+ ,x 2)") == List(:quasiquote,
                                  List(:(+), List(:unquote, :x), 2))
    @test_broken sx("`(+ ,@xs 2)") ==
        SExpression((:quasiquote,
                    (:(+),
                     (Symbol("unquote-splicing"), :xs),
                     2)))
end

@testset "file" begin
    @test SExpressions.parsefile("data/scheme.scm") isa List
end

@testset "macro" begin
    # Wrap in `@eval` to capture failures at test time rather than load time.
    @test @eval(sx"1")   == 1
    @test @eval(sx"foo") == :foo
end

@testset "interface" begin
    @test parse(SExpression, "(1 . + . (2 . * . 3))") == List(:+, 1, List(:*, 2, 3))
end

@testset "keywords" begin
    @test parse(SExpression, "(eq? '#:apple '#:orange)") ==
        List(Symbol("eq?"),
             List(:quote, Keyword("apple")),
             List(:quote, Keyword("orange")),
        )
end

end


text:
"""
    Basic <: TimeStrategy

Basic time strategy that adds new pedestrian to all available entrances every `at` iteration until `max_peds` pedestrians is added.

# Fields

- at::Int = 20
- max_peds::Int = 300
- vel_init::Point{<:Real} = initial velocity

"""
Base.@kwdef struct Basic{T} <: TimeStrategy
    at::Int = 20
    max_peds::Int = 300
    vel_init::Point{T} = (0.0, -1.3)
end

function add_pedestrian!(s::Basic, m::Model)
    if iter(m) == 0 || mod(iter(m), s.at) == 0
        for d in m.room.entrances
            maxid(m) >= s.max_peds && break
            add_pedestrian!(m, Pedestrian(0, middle(d); vel = s.vel_init))
        end
    end
    if maxid(m) == s.max_peds && npads(m) == 0
        set_stop!(m)
    end
    return
end


text:
struct DesignMatrix
    "Array of formulas"
    formulas::Any
    "A concatenated designmatric. In case of Mixed Models an array, where the first one is a FeMat, later ones ReMats. "
    Xs::Any
    "Event table with all events"
    events::Any
end

function DesignMatrix()
    return DesignMatrix([], [], [])
end

struct UnfoldMixedModelFitCollection{T<:AbstractFloat} <:
       MixedModels.MixedModelFitCollection{T}
    fits::Vector
    λ::Vector{<:Union{LowerTriangular{T,Matrix{T}},Diagonal{T,Vector{T}}}}
    inds::Vector{Vector{Int}}
    lowerbd::Vector{T}
    fcnames::NamedTuple
end

abstract type UnfoldModel end

#function UnfoldModel()
# todo make generator function with empty DesignMatrix
mutable struct UnfoldLinearModel <: UnfoldModel
    design::Dict
    designmatrix::DesignMatrix
    modelfit::Any
end

UnfoldLinearModel(d::Dict) = UnfoldLinearModel(d, Unfold.DesignMatrix(), [])
UnfoldLinearModel(d::Dict, X::DesignMatrix) = UnfoldLinearModel(d, X, [])

mutable struct UnfoldLinearMixedModel <: UnfoldModel
    design::Dict
    designmatrix::DesignMatrix
    modelfit::Any#::Array{UnfoldMixedModelFitCollection} # optional info on the modelfit
end
UnfoldLinearMixedModel(d::Dict) = UnfoldLinearMixedModel(d, Unfold.DesignMatrix(), [])
UnfoldLinearMixedModel(d::Dict, X::DesignMatrix) = UnfoldLinearMixedModel(d, X, [])

mutable struct UnfoldLinearModelContinuousTime <: UnfoldModel
    design::Dict
    designmatrix::DesignMatrix
    modelfit::Any
end

UnfoldLinearModelContinuousTime(d::Dict) =
    UnfoldLinearModelContinuousTime(d, Unfold.DesignMatrix(), [])
UnfoldLinearModelContinuousTime(d::Dict, X::DesignMatrix) =
    UnfoldLinearModelContinuousTime(d, X, [])

mutable struct UnfoldLinearMixedModelContinuousTime <: UnfoldModel
    design::Dict
    designmatrix::DesignMatrix
    modelfit::Any#::UnfoldMixedModelFitCollection
end

UnfoldLinearMixedModelContinuousTime(d::Dict) =
    UnfoldLinearMixedModelContinuousTime(d, Unfold.DesignMatrix(), [])
UnfoldLinearMixedModelContinuousTime(d::Dict, X::DesignMatrix) =
    UnfoldLinearMixedModelContinuousTime(d, X, [])

abstract type ModelFit end

struct LinearModelFit <: ModelFit
    estimate::Any
    info::Any
    standarderror::Any
end

LinearModelFit(estimate) = LinearModelFit(estimate, [], [])
LinearModelFit(estimate, info) = LinearModelFit(estimate, info, [])

function Base.show(io::IO, obj::UnfoldModel)
    println(io, "Unfold-Type: $(typeof(obj)) \n")
    println(io, "formula: $(obj.design)")
    println(
        io,
        "Useful functions:\n
    design(uf) \t\t(returns Dict of event => (formula,times/basis))  \n
    designmatrix(uf) \t(returns DesignMatrix with events) \n
    modelfit(uf) \t\t(returns modelfit object) \n
    coeftable(uf) \t\t(returns tidy result dataframe) \n",
    )
end

max_forks_repo_name: unfoldtoolbox/jlunfold


text:
import Base: show

mutable struct CompositePlot <: AbstractPlot
  plots::Vector{SimplePlot}
  size::Tuple

  function CompositePlot(plots, plot_size)
    @assert !is_repl "No composite plots in REPL"
    new(plots, plot_size)
  end
end

function CompositePlot(cur_plots=[]; kwargs...)
  cur_plot = nothing
  if haskey(kwargs, :figsize)
    cur_plot = CompositePlot(cur_plots, kwargs[:figsize])
  else
    cur_plot = CompositePlot(cur_plots, default_plot_size)
  end
  cur_plot
end


text:
###======================================================####
#Came from GLAbstraction/GLUniforms.jl

# Uniforms are OpenGL variables that stay the same for the entirety of a drawcall.
# There are a lot of functions, to upload them, as OpenGL doesn't rely on multiple dispatch.
# here is my approach, to handle all of the uniforms with one function, namely gluniform
# For uniforms, the Vector and Matrix types from ImmutableArrays should be used, as they map the relation almost 1:1

const GLSL_COMPATIBLE_NUMBER_TYPES = (GLfloat, GLint, GLuint, GLdouble)

function uniformfunc(typ::DataType, dims::Tuple{Int})
    Symbol(string("glUniform", first(dims), opengl_postfix(typ)))
end
function uniformfunc(typ::DataType, dims::Tuple{Int, Int})
    M, N = dims
    Symbol(string("glUniformMatrix", M == N ? "$M" : "$(M)x$(N)", opengl_postfix(typ)))
end

function gluniform(location::Integer, x::StaticArray)
    xref = [x]
    gluniform(location, xref)
end

@generated function gluniform(location::Integer, x::Vector{FSA}) where FSA <: StaticArray
    func = uniformfunc(eltype(FSA), size(FSA))
    callexpr = if ndims(FSA) == 2
        :($func(location, length(x), GL_FALSE, xref))
    else
        :($func(location, length(x), xref))
    end
    quote
        xref = reinterpret(eltype(FSA), x)
        $callexpr
    end
end

#Some additional uniform functions, not related to Imutable Arrays
# gluniform(location::Integer, target::Integer, t::Texture) = gluniform(GLint(location), GLint(target), t)
# gluniform(location::Integer, target::Integer, t::GPUVector) = gluniform(GLint(location), GLint(target), t.buffer)
# gluniform(location::Integer, target::Integer, t::TextureBuffer) = gluniform(GLint(location), GLint(target), t.texture)
# gluniform(location::Integer, t::TextureBuffer) = gluniform(GLint(location), GLint(target), t.texture)
#REVIEW: scary, binding and making texture active seems like something that shouldn't be in gluniform...

#TODO: We have the tools to figure out what texture unit a texture with a given name should be bound to, but we have to make it so that low level functionality remains to avoid slowness
gluniform(location::Integer, target::Integer, t::TextureBuffer) = gluniform(GLint(location), GLint(target), t.texture)
function gluniform(location::GLint, texture_unit, t::Texture)
    tu = GL_TEXTURE0 + UInt32(texture_unit)
    glActiveTexture(tu)
    glBindTexture(t.texturetype, t.id)
    gluniform(location, GLint(texture_unit))
end
gluniform(location::Integer, x::Enum)                                = gluniform(GLint(location), GLint(x))
gluniform(location::Integer, x::Union{GLubyte, GLushort, GLuint})    = glUniform1ui(GLint(location), x)
gluniform(location::Integer, x::Union{GLbyte, GLshort, GLint, Bool}) = glUniform1i(GLint(location),  x)
gluniform(location::Integer, x::GLfloat)                             = glUniform1f(GLint(location),  x)
gluniform(location::Integer, x::GLdouble)                            = glUniform1d(GLint(location),  x)

#Uniform upload functions for julia arrays...
gluniform(location::GLint, x::Vector{Float32})                       = glUniform1fv(location, length(x), x)
gluniform(location::GLint, x::Vector{GLdouble})                      = glUniform1dv(location, length(x), x)
gluniform(location::GLint, x::Vector{GLint})                         = glUniform1iv(location, length(x), x)
gluniform(location::GLint, x::Vector{GLuint})                        = glUniform1uiv(location, length(x), x)


text:
type MemcacheClients
    clients::Vector{MemcacheClient}

    function MemcacheClients(clients::Vector{MemcacheClient}=MemcacheClient[])
        new(clients)
    end
end

add_client(dmc::MemcacheClients, host::AbstractString, port::Integer=11211) = (push!(dmc.clients, MemcacheClient(host, port)); nothing)

for f in [:close, :quit, :flush_all, :slabs_automove, :slabs_reassign]
    @eval ($f)(dmc::MemcacheClients, args...) = (map(c->($f)(c, args...), dmc.clients); nothing)
end

for f in [:stats, :version]
    @eval ($f)(dmc::MemcacheClients, args...) = map(c->($f)(c, args...), dmc.clients)
end

mc(dmc::MemcacheClients, key::AbstractString) = dmc.clients[(hash(key) % length(dmc.clients)) + 1]
touch(dmc::MemcacheClients, key::AbstractString, exp::Integer; noreply::Bool=false) = touch(mc(dmc,key), key, exp, noreply=noreply)

for f in [:set, :add, :replace, :append, :prepend]
    @eval ($f)(dmc::MemcacheClients, key::AbstractString, val; exptime::Integer=0, cas::Integer=-1, noreply::Bool=false) = ($f)(mc(dmc,key), key, val; exptime=exptime, cas=cas, noreply=noreply)
end

cas(dmc::MemcacheClients, key::AbstractString, val, cas::Integer; exptime::Integer=0, noreply::Bool=false) = cas(mc(dmc,key), key, val, cas; exptime=exptime, noreply=noreply)

get(dmc::MemcacheClients, key::AbstractString) = get(mc(dmc,key), key)
get(dmc::MemcacheClients, key::AbstractString...; cas::Bool=false) = merge(map(c->get(c, key...; cas=cas), dmc.clients)...)
delete(dmc::MemcacheClients, key::AbstractString; noreply::Bool=false) = delete(mc(dmc,key), key, noreply=noreply)

incr(dmc::MemcacheClients, key::AbstractString, val::Integer; noreply::Bool=false) = incr(mc(dmc,key), key, val, noreply=noreply)
decr(dmc::MemcacheClients, key::AbstractString, val::Integer; noreply::Bool=false) = decr(mc(dmc,key), key, val, noreply=noreply)


text:
using RDatasets, KernelDensity
mlmf = dataset("mlmRev", "Gcsemv");

df = mlmf[complete_cases(mlmf[[:Written, :Course]]), :];
dc = array(df[:Course]);  kdc = kde(dc);
dw = array(df[:Written]); kdw = kde(dw);

using Winston
kdc = kde(dc); kdw = kde(dw);
plot(kdc.x, kdc.density, "r--", kdw.x, kdw.density, "b;")


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---
module SIRModel
using DifferentialEquations

function sir_ode(du, u, p, t)
    #Infected per-Capita Rate
    β = p[1]
    #Recover per-capita rate
    γ = p[2]

    #Susceptible Individuals
    S = u[1]
    #Infected by Infected Individuals
    I = u[2]

    du[1] = -β * S * I
    du[2] = β * S * I - γ * I
    du[3] = γ * I
end

#Pram = (Infected Per Capita Rate, Recover Per Capita Rate)
pram = [0.1,0.05]
#Initial Prams = (Susceptible Individuals, Infected by Infected Individuals)
init = [0.99,0.01,0.0]
tspan = (0.0,200.0)

sir_prob = ODEProblem(sir_ode, init, tspan, pram)

sir_sol = solve(sir_prob, saveat = 0.1);

# Visualization
using Plots
plot(sir_sol,xlabel="Time",ylabel="Number")

function sir_ode2(du,u,p,t)
    S,I,R = u
    b,g = p
    du[1] = -b*S*I
    du[2] = b*S*I-g*I
    du[3] = g*I
end
parms = [0.1,0.05]
init = [0.99,0.01,0.0]
tspan = (0.0,200.0)
sir_prob2 = ODEProblem(sir_ode2,init,tspan,parms)
sir_sol = solve(sir_prob2,saveat = 0.1)
end


text:
module Schedulers

using DataStructures

export
	Scheduler,
	JobStatusChange,
	PropagatedError,
	JobID,
	createjob!,
	deletejob!,
	add_dependency!,
	remove_dependency!,
	set_dependency!,
	setfunction!,
	setresult!,
	schedule!,
	cancel!,
	cancelall!,
	iscanceled,
	wantstorun,
	isactive,
	step!,
	status,
	statusstring,
	jobstatus,
	jobname

const JobID = Int
const Timestamp = Int
const Callback = Union{Function,Nothing}

struct JobStatus # TODO: rename
	changedAt::Threads.Atomic{Timestamp}
	runAt::Timestamp # i.e. the value of changedAt when the job was run
end
mutable struct Job
	changedAt::Threads.Atomic{Timestamp} # Only set by Scheduler thread. Can be read by other threads.
	name::String # Only used in error/log messages. Set once.
	f::Union{Function,Nothing}
	spawn::Bool
	edges::Dict{String,JobID} # name=>fromID
	edgesReverse::Set{Tuple{JobID,String}} # Set{(toID,name)}
	result::Any
	status::Symbol # one of :notstarted,:waiting,:spawned,:running,:done,:errored
	statusChangedTime::UInt64 # from time_ns()
	runAt::Timestamp
	waitingFor::Set{JobID}
	callbacks::Vector{Function}
end
Job(changedAt::Int, name::String, f::Union{Function,Nothing}, spawn::Bool, result, status::Symbol, runAt::Timestamp) = Job(Threads.Atomic{Int}(changedAt), name, f, spawn, Dict{String,JobID}(), Set{Tuple{JobID,String}}(), result, status, time_ns(), runAt, Set{JobID}(), [])

Job(name::String, f::Function, changedAt::Int, spawn::Bool) = Job(changedAt, name, f,       spawn, nothing, :notstarted, -1)
Job(name::String, result::Any, changedAt::Int, spawn::Bool) = Job(changedAt, name, nothing, spawn, result,  result2status(result), changedAt)

struct DetachedJob
	jobID::JobID
	runAt::Timestamp
end
struct JobStatusChange
	jobID::JobID
	jobName::String
	status::Symbol
	time::UInt64
	message::String
end
JobStatusChange(jobID::JobID,job::Job) =
	JobStatusChange(jobID, job.name, job.status, job.statusChangedTime, job.status==:errored ? sprint(showerror,job.result) : "")

struct Scheduler
	threaded::Bool
	timestamp::Ref{Timestamp}
	jobCounter::Threads.Atomic{JobID}
	jobs::Dict{JobID,Job}
	actions::Channel{Function}
	slowActions::Channel{Function}
	dirtyJobs::Vector{JobID}
	activeJobs::Set{JobID}
	detachedJobs::Dict{DetachedJob, UInt64} # -> time_ns() when job started running
	hasSpawned::Ref{Bool}
	statusChannel::Union{Channel{JobStatusChange},Nothing}
	verbose::Bool
end
function Scheduler(;threaded=Threads.nthreads()>1, statusChannel=nothing, verbose=false)
	threaded && Threads.nthreads()==1 && @warn "Trying to run threaded Scheduler, but threading is not enabled, please set the environment variable JULIA_NUM_THREADS to the desired number of threads."
	Scheduler(threaded, Ref(0), Threads.Atomic{JobID}(1), Dict{JobID,Job}(), Channel{Function}(Inf), Channel{Function}(Inf), [], Set{JobID}(), Dict{DetachedJob,UInt64}(), Ref(false), statusChannel, verbose)
end

struct PropagatedError{T<:Exception} <: Exception
	e::T
	jobName::String
	inputName::String
end
errorchain(io::IO, e::Exception) = showerror(io,e)
errorchain(io::IO, e::PropagatedError) = (print(io, e.inputName, '[', e.jobName, "]<--"); errorchain(io, e.e))
Base.showerror(io::IO, e::PropagatedError) = (print(io::IO, "Propagated error: "); errorchain(io,e))

# --- exported functions ---

function createjob!(s::Scheduler, result::Any, args::Pair{JobID,String}...; spawn::Bool=true, name::Union{String,Nothing}=nothing)
	id = newjobid!(s)
	name = name==nothing ? "AnonymousJob#$id" : name
	addaction!(s->_createjob!(s,id,result,name,spawn), s)
	for (fromID,toName) in args
		add_dependency!(s, fromID=>(id,toName))
	end
	id
end
createjob!(f::Function, s::Scheduler, args::Pair{JobID,String}...; name=string(f), kwargs...) = createjob!(s, f, args...; name=name, kwargs...)
deletejob!(s::Scheduler, jobID::JobID) = addaction!(s->_deletejob!(s,jobID), s)

add_dependency!(   s::Scheduler, dep::Pair{JobID,Tuple{JobID,String}}) = addaction!(s->   add_edge!(s,dep), s)
remove_dependency!(s::Scheduler, dep::Pair{JobID,Tuple{JobID,String}}) = addaction!(s->remove_edge!(s,dep), s)
set_dependency!(s::Scheduler, dep::Pair{JobID,Tuple{JobID,String}})    = addaction!(s->   set_edge!(s,dep), s)

setfunction!(f::Function, s::Scheduler, jobID::JobID) = addaction!(s->_setfunction!(s,jobID,f), s)
setresult!(s::Scheduler, jobID::JobID, result::Any)   = addaction!(s->_setresult!(s,jobID,result), s)

schedule!(callback::Callback, s::Scheduler, jobID::JobID) = addaction!(s->_schedule!(s,jobID,callback), s)
schedule!(s::Scheduler, jobID::JobID)                     = schedule!(nothing,s,jobID)

cancel!(s::Scheduler, jobID::JobID) = addaction!(s->_cancel!(s,jobID),s)
cancelall!(s::Scheduler) = addaction!(s) do s
	for jobID in keys(s.jobs)
		_cancel!(s, jobID)
	end
end

iscanceled(jobStatus::JobStatus) = jobStatus.changedAt[] > jobStatus.runAt

wantstorun(s::Scheduler) = isready(s.actions) || isready(s.slowActions)
function step!(s::Scheduler)
	while isready(s.actions)
		take!(s.actions)(s)
	end
	s.hasSpawned[] = false # TODO: revise this solution
	while !s.hasSpawned[] && isready(s.slowActions) # at most one slow spawning job per step
		take!(s.slowActions)(s)
	end
	updatetimestamps!(s)
end
isactive(s::Scheduler) = !isempty(s.activeJobs)

function status(s::Scheduler)
	# Running (A[2s], B[3s], ...), Detached Running (...), Spawned (C[0.1s], D [...), Waiting (...)
	now = time_ns()
	running  = Tuple{String,Float64}[]
	detached = Tuple{String,Float64}[]
	spawned  = Tuple{String,Float64}[]
	waiting  = Tuple{String,Float64}[]
	for jobID in s.activeJobs
		job = s.jobs[jobID]
		job.status == :running && push!(running,(job.name,(now-job.statusChangedTime)/1e9))
		job.status == :spawned && push!(spawned,(job.name,(now-job.statusChangedTime)/1e9))
		job.status == :waiting && push!(waiting,(job.name,(now-job.statusChangedTime)/1e9))
	end
	for (detachedJob,jobStartTime) in s.detachedJobs
		jobName = "Unknown"
		haskey(s.jobs, detachedJob.jobID) && (jobName = s.jobs[detachedJob.jobID].name)
		push!(detached, (jobName,(now-jobStartTime)/1e9))
	end
	running, detached, spawned, waiting
end
function statusstring(s::Scheduler; delim=", ")
	running, detached, spawned, waiting = status(s)
	strs = String[]
	isempty(running)  || push!(strs, string("Running (",  join(_durationstring.(running),  ", "), ")"))
	isempty(detached) || push!(strs, string("Detached (", join(_durationstring.(detached), ", "), ")"))
	isempty(spawned)  || push!(strs, string("Spawned (",  join(_durationstring.(spawned),  ", "), ")"))
	isempty(waiting)  || push!(strs, string("Waiting (",  join(_durationstring.(waiting),  ", "), ")"))
	join(strs, delim)
end

"""
	jobstatus(s::Scheduler, jobID::JobID)

Returns one of :doesntexist, :notstarted, :waiting, :spawned, :running, :done, :errored
"""
jobstatus(s::Scheduler, jobID::JobID) = haskey(s.jobs, jobID) ? s.jobs[jobID].status : :doesntexist

jobname(s::Scheduler, jobID::JobID) = haskey(s.jobs, jobID) ? s.jobs[jobID].name : "Unknown"

# --- internal functions ---

result2status(::Any)       = :done
result2status(::Exception) = :errored

function _setstatus!(s::Scheduler, jobID::JobID, job::Job, status::Symbol, statusChangedTime::UInt64=time_ns())
	if status in (:notstarted,:doesntexist) && job.status in (:running,:spawned)
		job.status==:running && @warn "Detaching running job $(job.name)"
		s.detachedJobs[DetachedJob(jobID,job.runAt)] = job.statusChangedTime
	end
	wasActive = job.status in (:waiting,:spawned,:running)
	isActive  =     status in (:waiting,:spawned,:running)
	job.status = status
	job.statusChangedTime = statusChangedTime
	wasActive && !isActive &&  pop!(s.activeJobs, jobID)
	isActive && !wasActive && push!(s.activeJobs, jobID)
	s.statusChannel!=nothing && put!(s.statusChannel, JobStatusChange(jobID, job))
	nothing
end
_durationstring(t::Tuple{String,Float64}, digits=1) = string(t[1], '[', round(t[2],digits=digits), "s]")

newjobid!(s::Scheduler) = Threads.atomic_add!(s.jobCounter, 1)
newtimestamp!(s::Scheduler) = s.timestamp[]+=1
addaction!(action::Function, s::Scheduler; slow::Bool=false) = (put!(slow ? s.slowActions : s.actions, action); nothing)

function _createjob!(s::Scheduler, jobID::JobID, x, name::String, spawn::Bool)
	job = Job(name, x, newtimestamp!(s), spawn)
	s.jobs[jobID] = job
	_setstatus!(s, jobID, job, job.status)
	nothing
end

function _deletejob!(s::Scheduler, jobID::JobID)
	@assert haskey(s.jobs, jobID) "Trying to delete nonexisting job with id $(jobID)"
	# remove dependencies
	job = s.jobs[jobID]
	for (name,fromID) in job.edges
		remove_edge!(s, fromID=>(jobID,name))
	end
	for (toID,name) in job.edgesReverse
		remove_edge!(s, jobID=>(toID,name))
	end
	_setstatus!(s, jobID, job, :doesntexist) # ensures removal from active jobs and detaching if needed
	delete!(s.jobs, jobID)
	nothing
end

function _setfunction!(s::Scheduler, jobID::JobID, f::Function)
	@assert haskey(s.jobs, jobID) "Trying to update nonexisting job with id $(jobID)"
	job = s.jobs[jobID]
	job.f = f
	setdirty!(s, jobID)
	nothing
end

function _setresult!(s::Scheduler, jobID::JobID, result::Any)
	@assert haskey(s.jobs, jobID) "Trying to update nonexisting job with id $(jobID)"
	job = s.jobs[jobID]
	@assert job.f == nothing
	result==job.result && return nothing # Nothing to do.
	setdirty!(s, jobID)
	updatetimestamps!(s)
	@assert job.status == :notstarted
	job.result = result
	_setstatus!(s, jobID, job, result2status(result))
	nothing
end

_schedule!(callback::Function, s::Scheduler, jobID::JobID, scheduledAt::Timestamp) = _schedule!(s,jobID,callback,scheduledAt)
function _schedule!(s::Scheduler, jobID::JobID, callback::Callback, scheduledAt::Timestamp = s.timestamp[])
	@assert haskey(s.jobs, jobID) "Trying to schedule nonexisting job with id $jobID"
	updatetimestamps!(s)
	job = s.jobs[jobID]
	changedAt = job.changedAt[]
	scheduledAt < changedAt && return # Scheduling out of date.
	slow = !job.spawn

	if job.runAt < changedAt # the job needs to run at a later timestamp than what has been run before
		@assert job.status == :notstarted

		job.runAt = changedAt
		_setstatus!(s, jobID, job, :waiting)
		empty!(job.callbacks)

		# figure out which jobs we are waiting for to finish
		waitingFor = Set{JobID}() # NB: do not reuse old set as it might be referenced from old dependency callbacks!
		for e in job.edges
			fromID = e[2]
			if s.jobs[fromID].result == nothing
				push!(waitingFor, fromID)
				_schedule!(s, fromID, scheduledAt) do x
					delete!(waitingFor,fromID)
					isempty(waitingFor) && addaction!(s->_spawn!(s,jobID,scheduledAt),s; slow=slow)
				end
			end
		end
		job.waitingFor = waitingFor
		s.verbose && @info "$(job.name) waiting for $(length(job.waitingFor)) jobs to finish."
	end
	@assert job.runAt == changedAt
	callback != nothing && push!(job.callbacks, callback)

	isempty(job.waitingFor) && addaction!(s->_spawn!(s,jobID,scheduledAt),s; slow=slow)

	nothing
end

function _setstatustorunning!(s::Scheduler, jobID::JobID, runAt::Timestamp, statusChangedTime::UInt64)
	haskey(s.jobs, jobID) || return # nothing to do if the job was deleted before it finished
	job = s.jobs[jobID]
	runAt < job.changedAt[] && return # nothing to do
	@assert job.status == :spawned
	_setstatus!(s, jobID, job, :running, statusChangedTime)
	nothing
end

function _spawn!(s::Scheduler, jobID::JobID, scheduledAt::Timestamp)
	haskey(s.jobs, jobID) || return # nothing to do if the job was deleted before it finished
	updatetimestamps!(s)
	job = s.jobs[jobID]
	changedAt = job.changedAt[]
	scheduledAt < changedAt && return # Spawned after dependency was finished, but now out of date.

	if job.status in (:done, :errored)
		job.result != nothing && _finish!(s, jobID, job.runAt, job.result, time_ns())
	elseif job.status == :waiting
		@assert isempty(job.waitingFor)
		s.verbose && @info "Spawning $(job.name) ($jobID@$changedAt)"
		s.hasSpawned[] = true
		inputs = getinput(s,jobID)
		_setstatus!(s, jobID, job, :spawned)
		f = job.f
		if s.threaded && job.spawn
			Threads.@spawn runjob!(s, jobID, job.name, job.changedAt, changedAt, f, inputs)
		else
			runjob!(s, jobID, job.name, job.changedAt, changedAt, f, inputs)
		end
	end
	nothing

end

function _cancel!(s::Scheduler, jobID::JobID)
	@assert haskey(s.jobs, jobID) "Trying to cancel nonexisting job with id $jobID"
	updatetimestamps!(s)
	job = s.jobs[jobID]
	job.result == nothing && setdirty!(s, jobID) # Do not invalidate jobs that have finished!
end

function _finish!(s::Scheduler, jobID::JobID, runAt::Timestamp, result::Any, statusChangedTime::UInt64)
	if haskey(s.jobs, jobID)
		updatetimestamps!(s)
		job = s.jobs[jobID]
		changedAt = job.changedAt[]
		if changedAt == runAt # we finished the last version of the job
			@assert job.status in (:running,:done,:errored)
			job.result = result
			_setstatus!(s, jobID, job, result2status(result), statusChangedTime)
			for callback in job.callbacks
				callback(job.result)
			end
			empty!(job.callbacks)
			result isa Exception && throw(result)
			return
		end
	end
	# if we reach here, the job was either deleted or we finished running an old version of the job
	pop!(s.detachedJobs, DetachedJob(jobID,runAt)) # remove from list of detached jobs
	nothing
end

function getinput(s::Scheduler, jobID::JobID)
	job = s.jobs[jobID]
	Dict{String,Any}((name=>s.jobs[inputJobID].result) for (name,inputJobID) in job.edges)
end

function add_edge!(s::Scheduler, dep::Pair{JobID,Tuple{JobID,String}})
	fromID,(toID,toName) = dep
	@assert haskey(s.jobs, toID)   "Trying to add edge to nonexisting job with id $(toID)"
	@assert haskey(s.jobs, fromID) "Trying to add edge from nonexisting job with id $(fromID)"
	from,to = s.jobs[fromID], s.jobs[toID]
	@assert !haskey(to.edges, toName) "Trying to add edge \"$toName\" that already exists in job with id ($toID)"
	to.edges[toName] = fromID
	push!(from.edgesReverse, (toID,toName))
	setdirty!(s, toID)
	nothing
end
function remove_edge!(s::Scheduler, dep::Pair{JobID,Tuple{JobID,String}})
	fromID,(toID,toName) = dep
	@assert haskey(s.jobs, toID)   "Trying to remove edge to nonexisting job with id $(toID)"
	@assert haskey(s.jobs, fromID) "Trying to remove edge from nonexisting job with id $(fromID)"
	from,to = s.jobs[fromID], s.jobs[toID]
	# we don't allow removing edges that doesn't exist
	@assert haskey(to.edges, toName) "Trying to remove edge \"$toName\" that doesn't exists in job with id ($toID)"
	@assert to.edges[toName]==fromID
	delete!(to.edges, toName)
	delete!(from.edgesReverse, (toID,toName))
	setdirty!(s, toID)
	nothing
end
function set_edge!(s::Scheduler, dep::Pair{JobID,Tuple{JobID,String}})
	fromID,(toID,toName) = dep
	@assert haskey(s.jobs, toID)   "Trying to replace edge to nonexisting job with id $(toID)"
	@assert haskey(s.jobs, fromID) "Trying to replace (remove) edge from nonexisting job with id $(fromID)"
	from,to = s.jobs[fromID], s.jobs[toID]
	if haskey(to.edges, toName)
		prevFromID = to.edges[toName]
		prevFromID == fromID && return
		remove_edge!(s, prevFromID=>(toID,toName))
	end
	add_edge!(s, dep)
	nothing
end

function setdirty!(s::Scheduler, jobID::JobID)
	job = s.jobs[jobID]
	job.changedAt[] = newtimestamp!(s)
	push!(s.dirtyJobs, jobID)
	nothing
end

function updatetimestamprec!(s::Scheduler, jobID::JobID)
	haskey(s.jobs, jobID) || return # happens if the job was deleted
	job = s.jobs[jobID]
	changedAt = job.changedAt[]

	_setstatus!(s, jobID, job, :notstarted)
	job.result = nothing

	for (toID,name) in job.edgesReverse
		toJob = s.jobs[toID]
		Threads.atomic_max!(toJob.changedAt,changedAt) < changedAt && updatetimestamprec!(s,toID)
	end
end

function updatetimestamps!(s::Scheduler)
	while !isempty(s.dirtyJobs) # go through in reverse order to minimize changes to dependency graph
		updatetimestamprec!(s,pop!(s.dirtyJobs))
	end
end

function runjob!(s::Scheduler, jobID::JobID, jobName::String, changedAt::Threads.Atomic{Timestamp}, runAt::Timestamp, f::Function, input::Dict{String,Any})
	jobStatus = JobStatus(changedAt, runAt)
	result = nothing
	if !iscanceled(jobStatus) # ensure the job was not changed after it was scheduled
		try
			s.verbose && @info "Running $jobName ($jobID@$runAt) with $(length(input)) parameter(s) in thread $(Threads.threadid())."
			jobStartTime = time_ns()
			addaction!(s->_setstatustorunning!(s,jobID,runAt,jobStartTime), s)

			# sanity check input and propagate errors
			for (name,v) in input
				@assert v != nothing "Job \"$jobName\" input \"$name\" has not been computed."
				v isa Exception && throw(PropagatedError(v,jobName,name))
			end
			result = f(jobStatus, input)
			dur = round((time_ns()-jobStartTime)/1e9,digits=1)
			isCanceled = iscanceled(jobStatus)
			detachedStr = isCanceled ? " (detached job)" : ""
			s.verbose && @info "Finished running$detachedStr $jobName[$(dur)s] ($jobID@$runAt) in thread $(Threads.threadid())."
			isCanceled || @assert result!=nothing "Job $jobName returned nothing"
		catch err
			if !(err isa PropagatedError)
				@warn "Job $jobName ($jobID@$runAt) failed"
				showerror(stdout, err, catch_backtrace())
			end
			result = err
		end
	end
	statusChangedTime = time_ns()
	addaction!(s->_finish!(s,jobID,runAt,result,statusChangedTime), s)
	result isa Exception && throw(result)
	nothing
end

end


text:
module TestUtils

using Test
using RLEVectors

@testset begin

    # rep
    @test RLEVectors.rep([4,5,6], each=2) == [4,4,5,5,6,6]
    @test RLEVectors.rep([4,5,6], times=2) == [4,5,6,4,5,6]
    @test RLEVectors.rep([4,5,6], each=3, times=2) == [4,4,4,5,5,5,6,6,6,4,4,4,5,5,5,6,6,6]
    @test RLEVectors.rep(["p","d","q"], each=3, times=2) == ["p","p","p","d","d","d","q","q","q","p","p","p","d","d","d","q","q","q"]
    @test_throws ArgumentError RLEVectors.rep([1,2,3], each=[4,5])
    @test rep(3,each=2) == [3,3]

    # searchsortedfirst
    @test RLEVectors.searchsortedfirst( [3, 6, 9, 12, 15 ], 6, 1, 5 ) == 2
    @test RLEVectors.searchsortedfirst( [3, 6, 9, 12, 15 ], 7, 1, 5 ) == 3
    @test RLEVectors.searchsortedfirst( [3, 6, 9, 12, 15 ], 7, 0, 22 ) == 3
    @test RLEVectors.searchsortedfirst( [3, 6, 9, 12, 15 ], 40, 1, 5 ) == 6
    @test RLEVectors.searchsortedfirst( [3, 6, 9, 12, 15 ], -2, 1, 5 ) == 1
    @test RLEVectors.searchsortedfirst( [3, 6, 9, 12, 15 ], 14, 4, 5 ) == 5
    @test RLEVectors.searchsortedfirst( [3, 6, 9, 12, 15 ], 3, 4, 5 ) == 4

    @test RLEVectors.searchsortedfirst( [0,5,10,15], [-3,2,3,7,22] ) == [1, 2, 2, 3, 5]
    @test RLEVectors.searchsortedfirst( [0,5,10,15], [-3,2,-3,7,22] ) == [1, 2, 1, 3, 5]

    # growat!
    x = collect(1:10)
    RLEVectors.growat!(x, 3, 2)
    @test x == [1,2,3,4,3,4,5,6,7,8,9,10]
    x = RLEVector([1,1,2,2,3,3])
    RLEVectors.growat!(x, 2, 1)
    @test x.runvalues[3:4] == [2,3]
    @test x.runends[3:4] == [4,6]

end # testset

end # module


text:
@testset "StructAxis" begin
    axis = @inferred(StructAxis{Complex{Float64}}())
    @test ArrayInterface.known_length(axis) === 2
    @test @inferred(getindex(axis, StaticInt(1):StaticInt(2))) === axis
    @test keys(@inferred(getindex(axis, StaticInt(1):StaticInt(1))))[1] === :re
    @test keys(@inferred(getindex(axis, StaticInt(2):StaticInt(2))))[1] === :im

    @testset "struct_view" begin
        x = AxisArray(reshape(1:4, 2, 2), StructAxis{Complex{Float64}}(), ["a", "b"]);
        @inferred(AxisIndices.structdim(x)) === StaticInt(1)
        @test_throws MethodError AxisIndices.structdim(parent(x))

        x2 = @inferred(struct_view(x))
        @test eltype(x2) <: Complex{Float64}
        @test keys(axes(x2, 1)) == ["a", "b"]

        x = AxisArray(reshape(1:4, 2, 2), StructAxis{NamedTuple{(:x, :y),Tuple{Float64,Float64}}}(), ["a", "b"]);
        x2 = @inferred(struct_view(x))
        @test eltype(x2) <: NamedTuple{(:x, :y),Tuple{Float64,Float64}}
        @test keys(axes(x2, 1)) == ["a", "b"]
    end

    # must know first index for struct axis b/c it's a completely statically sized axis
    @test_throws ArgumentError StructAxis{Complex{Float64}}(1:10)
    @test_throws ArgumentError StructAxis{Complex}(Base.OneTo(10))
end


text:
export PhasePartition
# Thermodynamic states
export ThermodynamicState,
       PhaseDry,
       PhaseEquil,
       PhaseNonEquil,
       TemperatureSHumEquil,
       LiquidIcePotTempSHumEquil,
       LiquidIcePotTempSHumNonEquil,
       LiquidIcePotTempSHumNonEquil_given_pressure

"""
    PhasePartition

Represents the mass fractions of the moist air mixture.

# Constructors

    PhasePartition(q_tot::Real[, q_liq::Real[, q_ice::Real]])
    PhasePartition(ts::ThermodynamicState)

See also [`PhasePartition_equil`](@ref)

# Fields

$(DocStringExtensions.FIELDS)
"""
struct PhasePartition{FT<:Real}
  "total specific humidity"
  tot::FT
  "liquid water specific humidity (default: `0`)"
  liq::FT
  "ice specific humidity (default: `0`)"
  ice::FT
end

PhasePartition(q_tot::FT,q_liq::FT) where {FT<:Real} =
  PhasePartition(q_tot, q_liq, zero(FT))
PhasePartition(q_tot::FT) where {FT<:Real} =
  PhasePartition(q_tot, zero(FT), zero(FT))

"""
    ThermodynamicState{FT}

A thermodynamic state, which can be initialized for
various thermodynamic formulations (via its sub-types).
All `ThermodynamicState`'s have access to functions to
compute all other thermodynamic properties.
"""
abstract type ThermodynamicState{FT} end

"""
    PhaseEquil{FT} <: ThermodynamicState

A thermodynamic state assuming thermodynamic equilibrium (therefore, saturation adjustment
may be needed).

# Constructors

    PhaseEquil(e_int, q_tot, ρ)

# Fields

$(DocStringExtensions.FIELDS)
"""
struct PhaseEquil{FT} <: ThermodynamicState{FT}
  "internal energy"
  e_int::FT
  "total specific humidity"
  q_tot::FT
  "density of air (potentially moist)"
  ρ::FT
  "temperature: computed via [`saturation_adjustment`](@ref)"
  T::FT
end
PhaseEquil(e_int::FT, q_tot::FT, ρ::FT) where {FT<:Real} =
  PhaseEquil(e_int, q_tot, ρ, saturation_adjustment(e_int, ρ, q_tot))

"""
    PhaseDry{FT} <: ThermodynamicState

A dry thermodynamic state (`q_tot = 0`).

# Constructors

    PhaseDry(e_int, ρ)

# Fields

$(DocStringExtensions.FIELDS)
"""
struct PhaseDry{FT} <: ThermodynamicState{FT}
  "internal energy"
  e_int::FT
  "density of dry air"
  ρ::FT
end

"""
    LiquidIcePotTempSHumEquil(θ_liq_ice, q_tot, ρ)

Constructs a [`PhaseEquil`](@ref) thermodynamic state from:

 - `θ_liq_ice` - liquid-ice potential temperature
 - `q_tot` - total specific humidity
 - `ρ` - density
"""
function LiquidIcePotTempSHumEquil(θ_liq_ice::FT, q_tot::FT, ρ::FT) where {FT<:Real}
    T = saturation_adjustment_q_tot_θ_liq_ice(θ_liq_ice, q_tot, ρ)
    q_pt = PhasePartition_equil(T, ρ, q_tot)
    e_int = internal_energy(T, q_pt)
    return PhaseEquil(e_int, q_tot, ρ, T)
end

"""
    TemperatureSHumEquil(T, q_tot, p)

Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.

 - `T` - temperature
 - `q_tot` - total specific humidity
 - `p` - pressure
"""
function TemperatureSHumEquil(T::FT, q_tot::FT, p::FT) where {FT<:Real}
    ρ = air_density(T, p, PhasePartition(q_tot))
    q = PhasePartition_equil(T, ρ, q_tot)
    e_int = internal_energy(T, q)
    PhaseEquil(e_int, q_tot, ρ, T)
end

"""
   	 PhaseNonEquil{FT} <: ThermodynamicState

A thermodynamic state assuming thermodynamic non-equilibrium (therefore, temperature can
be computed directly).

# Constructors

    PhaseNonEquil(e_int, q::PhasePartition, ρ)

# Fields

$(DocStringExtensions.FIELDS)

"""
struct PhaseNonEquil{FT} <: ThermodynamicState{FT}
  "internal energy"
  e_int::FT
  "phase partition"
  q::PhasePartition{FT}
  "density of air (potentially moist)"
  ρ::FT
end

"""
    LiquidIcePotTempSHumNonEquil(θ_liq_ice, q_pt, ρ)

Constructs a [`PhaseNonEquil`](@ref) thermodynamic state from:

 - `θ_liq_ice` - liquid-ice potential temperature
 - `q_pt` - phase partition
 - `ρ` - density
"""
function LiquidIcePotTempSHumNonEquil(θ_liq_ice::FT, q_pt::PhasePartition{FT}, ρ::FT) where {FT<:Real}
    T = air_temperature_from_liquid_ice_pottemp_non_linear(θ_liq_ice, ρ, q_pt)
    e_int = internal_energy(T, q_pt)
    PhaseNonEquil(e_int, q_pt, ρ)
end

"""
    LiquidIcePotTempSHumNonEquil_given_pressure(θ_liq_ice, q_pt, p)

Constructs a [`PhaseNonEquil`](@ref) thermodynamic state from:

 - `θ_liq_ice` - liquid-ice potential temperature
 - `q_pt` - phase partition
 - `p` - pressure
"""
function LiquidIcePotTempSHumNonEquil_given_pressure(θ_liq_ice::FT, q_pt::PhasePartition{FT}, p::FT) where {FT<:Real}
    T = air_temperature_from_liquid_ice_pottemp_given_pressure(θ_liq_ice, p, q_pt)
    ρ = air_density(T, p, q_pt)
    e_int = internal_energy(T, q_pt)
    PhaseNonEquil(e_int, q_pt, ρ)
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
mutable struct NetworkState <: State
    queues::Vector{Int} #A vector which indicates the number of customers in each queue
    in_transit::Int
    left_system::Int
    params::NetworkParameters #The parameters of the tandem queueing system
    arrival_distr::Gamma{Float64}
    travel_distr::Gamma{Float64}
    service_distrs::Vector{Gamma{Float64}}
end

"""
Constructs an initilized network state from the given parameters
"""
NetworkState(params::NetworkParameters) = NetworkState(
    fill(0, params.L), # Initial queues
    0,                  # Intial transit
    0,                  # Initial left
    params,
    rate_scv_gamma(params.λ, params.gamma_scv),
    rate_scv_gamma(params.η, params.gamma_scv),
    map((μ) -> rate_scv_gamma(μ, params.gamma_scv), params.μ_vector)
    )

queued_count(state::NetworkState)::Int = sum(state.queues)

transit_count(state::NetworkState)::Int = state.in_transit

num_in_queue(q::Int, state::NetworkState)::Int = state.queues[q]

function clear_left(state::NetworkState)
    state.left_system = 0
end

function show(io::IO, state::NetworkState)
    print(io, "State:\n")
    for (i, queue) in enumerate(state.queues)
        print(io, "Q$i: $queue\n")
    end
    print(io, "T: $(state.in_transit)\n")
end

function new_job(q::Int, time::Float64, state::NetworkState)
    state.queues[q] += 1
end

function new_transit(time::Float64, transit_time::Float64, state::NetworkState)
    state.in_transit += 1
end

function pop_queue(q::Int, time::Float64, state::NetworkState)
    state.queues[q] -= 1
    state.left_system += 1
end

function queue_to_transit(q::Int, time::Float64, transit_time::Float64, state::NetworkState)
    state.queues[q] -= 1
    state.in_transit += 1
end

function pop_transit(time::Float64, state::NetworkState)
    state.in_transit -= 1
    state.left_system += 1
end

function transit_to_queue(q::Int, time::Float64, state::NetworkState)
    state.in_transit -= 1
    state.queues[q] += 1
end

update_transit(time::Float64, transit_time::Float64, state::NetworkState) = nothing

function failed_arrival(time::Float64, state::NetworkState)
    state.left_system += 1
end


text:
export SnakeGame

using Random

struct SnakeGame{N,M,R<:AbstractRNG}
    # WHCN
    board::BitArray{M}
    # TODO: using DataStructures: Deque ?
    snakes::Vector{Vector{CartesianIndex{N}}}
    foods::Set{CartesianIndex{N}}
    walls::Vector{CartesianIndex{N}}
    n_foods::Int
    rng::R
end

Base.size(g::SnakeGame) = size(g.board)[1:end-1]

ind_of_wall(game) = 2*length(game.snakes)+1
ind_of_food(game) = 2*length(game.snakes)+2

mark_wall!(game, ind) = game.board[ind, ind_of_wall(game)] = true
mark_food!(game, ind) = game.board[ind, ind_of_food(game)] = true
unmark_food!(game, ind) = game.board[ind, ind_of_food(game)] = false
mark_snake_head!(game, ind, i) = game.board[ind, i] = true
unmark_snake_head!(game, ind, i) = game.board[ind, i] = false
mark_snake_body!(game, ind, i) = game.board[ind, length(game.snakes)+i] = true
unmark_snake_body!(game, ind, i) = game.board[ind, length(game.snakes)+i] = false

get_walls(game) = findall(isone, selectdim(game.board, ndims(game.board), ind_of_wall(game)))

"""
    SnakeGame(;kwargs...)

# Keyword Arguments

- `shape::NTuple{N,Int}=(8,8)`, the size of game board. `N` can be greater than 2.
- `walls`, an iterable type with elements of type `CartesianIndex{N}`.
- `n_snakes`, number of snakes.
- `n_foods`, maximum number of foods in each step.
- `rng::AbstractRNG`, inner RNG used to sample initial snakes and necessary foods in each step.
"""
function SnakeGame(;shape=(8,8), walls=CartesianIndex{length(shape)}[], n_snakes=1, n_foods=1,rng=Random.GLOBAL_RNG)
    n_snakes+n_foods >= reduce(*, shape) && error("n_snakes+n_foods must be less than the total grids")
    board = BitArray(undef, shape..., n_snakes#=snake head=#+n_snakes#=snake body=#+1#=wall=#+1#=food=#)
    snakes = [Vector{CartesianIndex{length(shape)}}() for _ in 1:n_snakes]
    foods = Set{CartesianIndex{length(shape)}}()
    game = SnakeGame(board, snakes, foods, walls, n_foods, rng)
    reset!(game)
    game
end

function reset!(game::SnakeGame)
    fill!(game.board, false)

    for s in game.snakes
        empty!(s)
    end

    empty!(game.foods)

    # do not change walls

    for w in game.walls
        mark_wall!(game, w)
    end

    while length(game.foods) < game.n_foods
        p = rand(game.rng, CartesianIndices(size(game)))
        if any(@view(game.board[p, :]))
            continue  # there's a wall
        else
            push!(game.foods, p)
            mark_food!(game, p)
        end
    end

    for i in 1:length(game.snakes)
        while true
            p = rand(game.rng, CartesianIndices(size(game)))
            if any(@view(game.board[p, :]))
                continue  # there's a wall or food
            else
                push!(game.snakes[i], p)
                mark_snake_head!(game, p, i)
                break
            end
        end
    end

end

(game::SnakeGame)(action::CartesianIndex) = game([action])

function (game::SnakeGame{N})(actions::Vector{CartesianIndex{N}}) where N
    # 1. move snake
    for ((i, s), a) in zip(enumerate(game.snakes), actions)
        unmark_snake_head!(game, s[1], i)
        mark_snake_body!(game, s[1], i)
        pushfirst!(s, CartesianIndex(mod.((s[1] + a).I, axes(game.board)[1:end-1])))
        mark_snake_head!(game, s[1], i)
        if s[1] in game.foods
            unmark_food!(game, s[1])
        else
            unmark_snake_body!(game, s[end], i)
            pop!(s)
        end
    end
    # 2. check collision
    for s in game.snakes
        sum(@view(game.board[s[1],:])) == 1 || return false
    end
    # 3. create new foods
    for f in game.foods
        if !game.board[f, ind_of_food(game)]
            # food is eaten
            pop!(game.foods, f)
            food = rand(game.rng, CartesianIndices(size(game)))
            attempts = 1
            while any(@view(game.board[food, :]))
                food = rand(game.rng, CartesianIndices(size(game)))
                attempts += 1
                if attempts > reduce(*, size(game))
                    @warn "a rare case happened: sampled too many times to generate food"
                    empty_positions = findall(iszero, vec(any(game.board, dims=ndims(game.board))))
                    if length(empty_positions) == 0
                        return false
                    else
                        food = CartesianIndices(size(game))[rand(game.rng, empty_positions)]
                        break
                    end
                end
            end
            push!(game.foods, food)
            mark_food!(game, food)
        end
    end
    return true
end


text:
# 00-data.jl
# Extracts audio features from TIMIT to be used in speech recognition

using Flux: onehotbatch
using WAV
using BSON

# This wookay's fork of MFCC updated to work with Julia v0.7/1.0
# https://github.com/wookay/MFCC.jl
using MFCC

# Define constants that will be used
const TRAINING_DATA_DIR = "TIMIT/TRAIN"
const TEST_DATA_DIR = "TIMIT/TEST"

const TRAINING_OUT_DIR = "train"
const TEST_OUT_DIR = "test"

# Make dictionary to map from phones to class numbers
const PHONES = split("h#	q	eh	dx	iy	r	ey	ix	tcl	sh	ow	z	s	hh	aw	m	t	er	l	w	aa	hv	ae	dcl	y	axr	d	kcl	k	ux	ng	gcl	g	ao	epi	ih	p	ay	v	n	f	jh	ax	en	oy	dh	pcl	ah	bcl	el	zh	uw	pau	b	uh	th	ax-h	em	ch	nx	eng")
translations = Dict(phone=>i for (i, phone) in enumerate(PHONES))
translations["sil"] = translations["h#"]
const PHONE_TRANSLATIONS = translations

# Make dictionary to perform class folding
const FOLDINGS = Dict(
  "ao" => "aa",
  "ax" => "ah",
  "ax-h" => "ah",
  "axr" => "er",
  "hv" => "hh",
  "ix" => "ih",
  "el" => "l",
  "em" => "m",
  "en" => "n",
  "nx" => "n",
  "eng" => "ng",
  "zh" => "sh",
  "pcl" => "sil",
  "tcl" => "sil",
  "kcl" => "sil",
  "bcl" => "sil",
  "dcl" => "sil",
  "gcl" => "sil",
  "h#" => "sil",
  "pau" => "sil",
  "epi" => "sil",
  "ux" => "uw"
)

FRAME_LENGTH = 0.025 # ms
FRAME_INTERVAL = 0.010 # ms

"""
  makeFeatures(wavFname, phnFname)

Extracts Mel filterbanks and associated labels from `wavFname` and `phnFaname`.
"""
function makeFeatures(phnFname, wavFname)
  samps, sr = wavread(wavFname)
  samps = vec(samps)

  mfccs, _, _ = mfcc(samps, sr, :rasta; wintime=FRAME_LENGTH, steptime=FRAME_INTERVAL)

  local lines
  open(phnFname, "r") do f
    lines = readlines(f)
  end

  boundaries = Vector()
  labels = Vector()

  # first field in the file is the beginning sample number, which isn't
  # needed for calculating where the labels are
  for line in lines
    _, boundary, label = split(line)
    boundary = parse(Int64, boundary)
    push!(boundaries, boundary)
    push!(labels, label)
  end

  labelInfo = collect(zip(boundaries, labels))
  labelInfoIdx = 1
  boundary, label = labelInfo[labelInfoIdx]
  nSegments = length(labelInfo)

  frameLengthSamples = FRAME_LENGTH * sr
  frameIntervalSamples = FRAME_INTERVAL * sr
  halfFrameLength = FRAME_LENGTH / 2

  # Begin generating sequence labels by looping through the MFCC
  # frames

  labelSequence = Vector() # Holds the sequence of labels

  idxsToDelete = Vector() # To store indices for frames labeled as 'q'
  for i=1:size(mfccs, 1)
    win_end = frameLengthSamples + (i-1)*frameIntervalSamples

    # Move on to next label if current frame of samples is more than half
    # way into next labeled section and there are still more labels to
    # iterate through
    if labelInfoIdx < nSegments && win_end - boundary > halfFrameLength

      labelInfoIdx += 1
      boundary, label = labelInfo[labelInfoIdx]
    end

    if label == "q"
      push!(idxsToDelete, i)
      continue
    end

    push!(labelSequence, label)
  end

  # Remove the frames that were labeld as 'q'
  mfccs = mfccs[[i for i in 1:size(mfccs,1) if !(i in Set(idxsToDelete))],:]

  mfccDeltas = deltas(mfccs, 2)
  features = hcat(mfccs, mfccDeltas)
  return (features, labelSequence)
end

"""
  createData(data_dir, out_dir)

Extracts data from files in `data_dir` and saves results in `out_dir`.
"""
function createData(data_dir, out_dir)

  ! isdir(out_dir) && mkdir(out_dir)

  for (root, dirs, files) in walkdir(data_dir)

    # Exclude the files that are part of the speaker accent readings
    files = [x for x in files if ! occursin("SA", x)]

    phnFnames = [x for x in files if occursin("PHN", x)]
    wavFnames = [x for x in files if occursin("WAV", x)]

    one_dir_up = basename(root)
    print("$(root)\r")

    for (wavFname, phnFname) in zip(wavFnames, phnFnames)
      phn_path = joinpath(root, phnFname)
      wav_path = joinpath(root, wavFname)

      x, y = makeFeatures(phn_path, wav_path)

      # Generate class nums; there are 61 total classes, but only 39 are
      # used after folding.
      y = [PHONE_TRANSLATIONS[x] for x in y]
      class_nums = [n for n in 1:61]
      y = onehotbatch(y, class_nums)

      base, _ = splitext(phnFname)
      dat_name = one_dir_up * base * ".bson"
      dat_path = joinpath(out_dir, dat_name)
      BSON.@save dat_path x y
    end
  end
  println()
end

createData(TRAINING_DATA_DIR, TRAINING_OUT_DIR)
createData(TEST_DATA_DIR, TEST_OUT_DIR)

max_stars_repo_name: jstrube/model-zoo

max_issues_repo_name: jstrube/model-zoo

max_forks_repo_name: jstrube/model-zoo


text:
# ======================
# Dependencies
# ======================

using LinearAlgebra, SparseArrays, # modules from the Julia standard library
      Reexport,                    # see @reexport macro below
      RecipesBase,                 # plotting
      Parameters,                  # structs with kwargs
      StaticArrays,                # statically sized arrays
      RecursiveArrayTools          # vector of arrays type

# manipulate function definition expressions
using ExprTools: splitdef, combinedef

# the reexport macro ensures that the names exported by the following libraries
# are made available after loading ReachabilityAnalysis
@reexport using HybridSystems,
                IntervalMatrices,
                LazySets,
                MathematicalSystems,
                TaylorIntegration

# required to avoid conflicts with MathematicalSystems
using LazySets: LinearMap, AffineMap, ResetMap

# required to avoid conflicts with IntervalMatrices
using LazySets: Interval, radius, sample, ∅, dim, scale, scale!

# in-place set operations
using LazySets: linear_map!

# LazySets internal functions frequently used
using LazySets.Arrays: projection_matrix, SingleEntryVector
using LazySets.Approximations: AbstractDirections
using LazySets: @commutative

# aliases for intervals
const IM = IntervalMatrices
import IntervalArithmetic
const IA = IntervalArithmetic
const TimeInterval = IA.Interval{Float64}
import TaylorModels
const TM = TaylorModels

# method extensions for Taylor model reach-sets
import TaylorModels: domain, remainder, polynomial, get_order, evaluate

# aliases for set types
const CPA = CartesianProductArray

# convenience union for dispatch on structs that are admissible as initial sets or inputs
const AdmissibleSet = Union{LazySet, UnionSet, UnionSetArray, IA.Interval, IA.IntervalBox}

# method extensions
import LazySets: dim, overapproximate, box_approximation, project, Projection,
                 intersection, directions, linear_map, LinearMap, _split, split!,
                 set, array, _isapprox

import MathematicalSystems: discretize
import Base: ∈, ∩, convert, isdisjoint
import LinearAlgebra: normalize

import CommonSolve: solve # common solve name

# ======================
# Useful constants
# ======================

@inline zeroBox(m) = IntervalBox(zeroI, m)
@inline unitBox(m) = IntervalBox(IA.Interval(0.0, 1.0), m)
@inline symBox(n::Integer) = IntervalBox(symI, n)
const zeroI = IA.Interval(0.0) # TODO use number type
const oneI = IA.Interval(1.0)
const symI = IA.Interval(-1.0, 1.0)

# Interval constructor given a float
LazySets.Interval(x::Float64) = Interval(interval(x))

# common aliases for system's names
const LCS = LinearContinuousSystem
const LDS = LinearDiscreteSystem
const CLCS = ConstrainedLinearContinuousSystem
const CLDS = ConstrainedLinearDiscreteSystem
const CLCCS = ConstrainedLinearControlContinuousSystem
const CLCDS = ConstrainedLinearControlDiscreteSystem
const ACS = AffineContinuousSystem
const ADS = AffineDiscreteSystem
const CACCS = ConstrainedAffineControlContinuousSystem
const CACDS = ConstrainedAffineControlDiscreteSystem
const CACS = ConstrainedAffineContinuousSystem
const CADS = ConstrainedAffineDiscreteSystem
const BBCS = BlackBoxContinuousSystem
const CBBCS = ConstrainedBlackBoxContinuousSystem
const CBBCCS = ConstrainedBlackBoxControlContinuousSystem
const SOLCS = SecondOrderLinearContinuousSystem
const SOACS = SecondOrderAffineContinuousSystem
const SOCLCCS = SecondOrderConstrainedLinearControlContinuousSystem
const SOCACCS = SecondOrderConstrainedAffineControlContinuousSystem
const SecondOrderSystem = Union{SOLCS, SOACS, SOCLCCS, SOCACCS}
const NonlinearSystem = Union{BBCS, CBBCS, CBBCCS}

@inline function _isapprox(Δt::TimeInterval, Δs::TimeInterval)
    return (inf(Δt) ≈ inf(Δs)) && (sup(Δt) ≈ sup(Δs))
end

const VecOrTuple = Union{<:AbstractVector{Int}, NTuple{D, Int}} where {D}
const VecOrTupleOrInt = Union{<:AbstractVector{Int}, NTuple{D, Int}, Int} where {D}

# ======================
# Optional dependencies
# ======================

using Requires

function __init__()
    # numerical differential equations suite
    @require DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" include("init_DifferentialEquations.jl")

    # exponentiation methods using Krylov subspace approximations
    @require ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" include("init_ExponentialUtilities.jl")

    # methods that require optimization modeling and solvers
    @require JuMP = "4076af6c-e467-56ae-b986-b466b2749572" include("init_JuMP.jl")

    # tools for symbolic computation
    #@require ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" include("init_ModelingToolkit.jl")
    @require Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" include("init_Symbolics.jl")

    # tools for symbolic algebra
    @require MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" include("init_MultivariatePolynomials.jl")

    # sparse dynamic representation of multivariate polynomials
    @require DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" include("init_DynamicPolynomials.jl")
end

# ===========================
# Utility macros
# ===========================

"""
    @requires(module_name)

Convenience macro to annotate that a package is required to use a certain function.

### Input

- `module_name` -- name of the required package

### Output

The macro expands to an assertion that checks whether the module `module_name` is
known in the calling scope.

### Notes

Usage:

```julia
function foo(...)
    @require MyPackage
    ... # functionality that requires MyPackage to be loaded
end
```
"""
macro requires(module_name)
    m = Meta.quot(Symbol(module_name))
    return esc(:(@assert isdefined(@__MODULE__, $m) "package `$($m)` is required " *
                    "for this function; do `using $($m)` and try again"))
end


text:
mutable struct Translational{T,Nc} <: Joint{T,Nc}
    vertices::NTuple{2,SVector{3,T}}
    V12::SMatrix{2,3,T,6}
    V3::Adjoint{T,SVector{3,T}}
    cid::Int64

    function Translational{T,Nc}(body1::AbstractBody{T}, body2::AbstractBody{T};
        p1::AbstractVector{T} = zeros(3), p2::AbstractVector{T} = zeros(3), axis::AbstractVector{T} = zeros(3)) where {T,Nc}

        vertices = (p1, p2)
        if norm(axis) != 0
            axis = axis / norm(axis)
        end
        A = svd(skew(axis)).Vt # in frame of body1
        V12 = A[1:2,:]
        V3 = axis' # instead of A[3,:] for correct sign: abs(axis) = abs(A[3,:])
        cid = body2.id

        new{T,Nc}(vertices, V12, V3, cid), body1.id, body2.id
    end
end

Translational0 = Translational{T,0} where T
Translational1 = Translational{T,1} where T
Translational2 = Translational{T,2} where T
Translational3 = Translational{T,3} where T


text:
module RayTracingWeekend

using Images
using LinearAlgebra
using Random
using RandomNumbers.Xorshifts
using StaticArrays

include("vec.jl")
export Vec3, squared_length, near_zero, rgb, rgb_gamma2
include("init.jl")
export TRNG
include("structs.jl")
export Ray, Hittable, HittableList, Material, HitRecord, Sphere, Scatter
include("rand.jl")
export reseed!, trand, random_vec3_in_sphere, random_between, random_vec3, random_vec2,
	random_vec3_on_sphere, random_vec2_in_disk
include("hit.jl")
export point, ray_to_HitRecord, hit
include("ray_color.jl")
export skycolor, color_vec3_in_rgb, ray_color
include("camera.jl")
export Camera, default_camera, get_ray
include("render.jl")
export render
include("light.jl") # light transport
export reflect, refract, reflectance
include("material.jl")
export Lambertian, scatter, Metal, Dielectric
include("scenes.jl")
export scene_2_spheres, scene_4_spheres, scene_blue_red_spheres, scene_diel_spheres, scene_random_spheres

end


text:
#=
Root file for the routing related files.
=#

module Routing

using ..Mapper2.Helper
Helper.@SetupDocStringTemplates

using ..Mapper2.MapperCore
using ..Mapper2.MapperGraphs

using DataStructures

using Logging
using Random

export route!,
    ChannelIndex,
    PortVertices,
    RoutingStruct,
    RoutingLink,
    BasicRoutingLink,
    RoutingChannel,
    BasicRoutingChannel,
    routing_link_type,
    routing_channel,
    annotate

abstract type AbstractRoutingAlgorithm end

"""
Representation of routing resources in an architecture.

API
---
* [`channels`](@ref)
* [`cost`](@ref)
* [`capacity`](@ref)
* [`occupancy`](@ref)
* [`addchannel`](@ref)
* [`remchannel`](@ref)

Implementations
---------------
* [`BasicRoutingLink`](@ref) - Reference this type for what methods of the API
    come for free when using various fields of the basic type.
"""
abstract type RoutingLink end

"""
Representation of channels in the taskgraph for routing.

API
---
* [`start_vertices`](@ref)
* [`stop_vertices`](@ref)
* [`isless(a::RoutingChannel, b::RoutingChannel)`](@ref)

Implementations
---------------
* [`BasicChannel`](@ref)
"""
abstract type RoutingChannel end

################################################################################
include("Types.jl")
include("Graph.jl")
include("Links.jl")
include("Channels.jl")
include("Struct.jl")

# Algorithms
include("Pathfinder/Pathfinder.jl")

# Routing dispatch
"""
    route!(map::Map; meta_prefix = "")

Run pathfinder routing directly on `map`. Keyword argument `meta_prefix` allows
controlling the prefix of the `map.metadata` keys generated below.

Records the following metrics into `map.metadata`:

* `routing_struct_time` - Time it took to build the [`RoutingStruct`](@ref)

* `routing_struct_bytes` - Memory allocation to build [`RoutingStruct`](@ref)

* `routing_passed :: Bool` - `true` if routing passes [`check_routing`](@ref),
    otherwise `false`.

* `routing_error :: Bool` - `true` if routing experienced a connectivity error.

The following are also included if `routing_error == false`.

* `routing_time` - Time to run routing.

* `routing_bytes` - Memory allocation during routing.

* `routing_global_links` - The number of global links used in the final routing.
"""
function route!(m::Map; meta_prefix = "")
    # Build the routing structure
    routing_struct, struct_time, struct_bytes, _, _ = @timed RoutingStruct(m)
    # Default to Pathfinder
    algorithm = routing_algorithm(m, routing_struct)
    # Run the routing algorithm.
    routing_error = false
    local route_time
    local route_bytes
    try
        _, route_time, route_bytes, _, _ = @timed route!(algorithm, routing_struct)
    catch err
        @error err
        routing_error = true
    end

    # Record the final results.
    record(m, routing_struct)
    routing_passed = check_routing(m)

    # Save all of this to metadata.
    m.metadata["$(meta_prefix)routing_struct_time"] = struct_time
    m.metadata["$(meta_prefix)routing_struct_bytes"] = struct_bytes
    m.metadata["$(meta_prefix)routing_passed"] = routing_passed
    m.metadata["$(meta_prefix)routing_error"] = routing_error
    if !routing_error
        m.metadata["$(meta_prefix)routing_time"] = route_time
        m.metadata["$(meta_prefix)routing_bytes"] = route_bytes
        m.metadata["$(meta_prefix)routing_global_links"] = MapperCore.total_global_links(m)
    end

    return nothing
end

routing_algorithm(m::Map, rs) = Pathfinder(m, rs)

################################################################################
# REQUIRED METHODS
################################################################################
"""
    annotate(ruleset::RuleSet, item::Union{Port,Link,Component})

Return some [`<:RoutingLink`](@ref RoutingLink) for `item`.  If
`item <: Component`, it is a primitive. If not other primitives have been
defined, it will be a `mux`.

See: [`BasicRoutingLink`](@ref)

Default: `BasicRoutingLink(capacity = getcapacity(ruleset, item))
"""
function annotate(ruleset::RuleSet, item::Union{Port,Link,Component})
    return BasicRoutingLink(; capacity = getcapacity(ruleset, item))
end

"""
    canuse(ruleset::RuleSet, link::RoutingLink, channel::RoutingChannel)::Bool

Return `true` if `channel` can be routed using `link`.

See: [`RoutingLink`](@ref), [`RoutingChannel`](@ref)

Default: `true`
"""
MapperCore.canuse(::RuleSet, link::RoutingLink, channel::RoutingChannel) = true

end


text:
# Should find more elgant method to find required atype
_get_basetype(x::Type{Knet.KnetArray{T,N}}) where {T,N} = Knet.KnetArray{T}
_get_basetype(x::Type{Knet.AutoGrad.Result{Knet.KnetArray{T,N}}}) where {T,N} = Knet.KnetArray{T}
_get_basetype(x::Type{Array{T,N}}) where {T,N} = Array{T}
_get_basetype(x::Type{Knet.AutoGrad.Result{Array{T,N}}}) where {T,N} = Array{T}

function getparam(atype::Type, d...;init=Knet.xavier)
	atype = _get_basetype(atype)
	Knet.Param(atype(init(d...)))
end

const RegisteredActivations = IdDict{String, Any}(
	"relu" => Knet.relu,
	"elu" => Knet.elu,
	"selu" => Knet.selu,
	"sigm" => Knet.sigm,
	"sigmoid" => Knet.sigm,
	"identity" => identity
)

get_activation(f) = f
function get_activation(f::Union{String, Symbol})
	f = string(f)
	if ! haskey(RegisteredActivations, f)
		@warn "unkown activtion function, using relu instead" fn=f
		return Knet.relu
	else
		return RegisteredActivations[f]
	end
end

const _layernames = Dict{String,Int}()

function get_layername(layername; kwargs...)
	info = Dict(kwargs)
	if haskey(info, :name)
		info[:name]
	else
		if ! haskey(_layernames, layername)
			_layernames[layername] = 1
		end
		id = _layernames[layername]
		_layernames[layername] += 1
		"$(layername)_$(id)"
	end
end

abstract type Layer end

"""
Lazy layers will be created during the first invocation. This means that a development time
you don't need to specify much, since most aspects will be figured out based on the past data.
"""
abstract type LazyLayer <: Layer end

## Generic code, X can be a Tensor or a Tuple
function (layer::Layer)(X)
	@debug "Calling $(typeof(layer))" layer
    call(layer, X)
end

## Generic code, X can be a Tensor or a Tuple
function (layer::LazyLayer)(X)
	@assert hasproperty(layer, :built) "LazyLayer $layer without built property detected"
	@assert isa(layer.built, Bool) "LazyLayer $layer built property not a boolean"

    if ! layer.built
		@debug "Building $(typeof(layer))" layer
        build(layer, size(X)[1:end-1], typeof(X))
		layer.built = true
    end
	@debug "Calling $(typeof(layer))" layer
    call(layer, X)
end

"""
Regular densely-connected NN layer.

Dense implements the operation: output = activation(dot(input, weight) + bias) where
activation is the element-wise activation function passed as the activation argument,
weight is a weights matrix created by the layer, and bias is a bias vector created
by the layer (only applicable if use_bias is true).

# Usage

```julia
layer = Dense(10, relu)
layer = Dense(100, use_bias=false)
```

"""
mutable struct Dense <:LazyLayer
    units::Int
	activation::Function
	use_bias::Bool
	name::String
	built::Bool
	init::NamedTuple
	params::NamedTuple

	function Dense(units::Int, activation=identity; use_bias=true,
		initw = Knet.xavier, initb = zeros, kwargs...)
		@assert units > 0 "Units of a Dense layer should be > 0"
		activation = get_activation(activation)
		name = get_layername("dense"; kwargs...)
	    new(units, activation, use_bias, name, false,
		(w=initw, b=initb), (w=nothing, b=nothing))
	end
end

function call(layer::Dense, X::Tensor)
	X = Knet.mat(X) # Flatten if required

	w,b = layer.params
	if layer.use_bias
    	layer.activation.(w*X .+ b)
	else
		layer.activation.(w*X)
	end
end

function build(layer::Dense, shape::Shape, atype)
	# nInput = length(shape) > 1 ? *(shape...) : shape[1]
	nInput = *(shape...)

    w = getparam(atype, layer.units, nInput, init=layer.init.w)
	b = nothing
    if layer.use_bias
		b = getparam(atype, layer.units, init=layer.init.b)
	end
	layer.params = (w=w,b=b)
end

"""
Softmax
"""
mutable struct Softmax <: Layer
	dims
	Softmax(dims=:) = new(dims)
end

function call(layer::Softmax, X::Tensor)
	Knet.softmax(X, dims=layer.dims)
end

"""
Log Softmax
"""
mutable struct LogSoftmax <: Layer
	dims
	LogSoftmax(dims=:) = new(dims)
end

function call(layer::LogSoftmax, X::Tensor)
	Knet.logsoftmax(X, dims=layer.dims)
end

"""
Flattening Layer. Photon by default already has flattening functionality
build into the Dense layer, so you won't need to include a separate Flatten
layer before a Dense layer.
"""
mutable struct Flatten <: Layer
	dims
	Flatten(dims=nothing) = new(dims)
end

function call(layer::Flatten, X::Tensor)
	layer.dims === nothing ? Knet.mat(X) : Knet.mat(X, dims=layer.dims)
end

"""
Reshape the shape of a tensor.
"""
mutable struct Reshape <: Layer
	dims
	Reshape(dims::Tuple) = new(dims)
end

function call(layer::Reshape, X::Tensor)
	dims = collect(layer.dims)
    batchsize = size(X)[end]
    newdim = push!(dims, batchsize)
	reshape(X, Tuple(newdim))
end

## Activation
mutable struct Activation <:Layer
	activation
end

function call(c::Activation, X::Tensor)
	c.activation.(X)
end

"""
Dropout layer with optional the rate (between 0 and 1) of dropout. If
no rate is specified, 0.5 (so 50%) will be used.
"""
struct Dropout <: Layer
	rate::Float64

	function Dropout(rate=0.5)
		@assert rate <= 1.0
		@assert rate >= 0.0
		new(rate)
	end
end

function call(d::Dropout, X::Tensor)
	Knet.dropout(X, d.rate)
end

"""
Batch normalization layer (Ioffe and Szegedy, 2014). Normalizes the input at each batch,
i.e. applies a transformation that maintains the mean activation close to 0
and the activation standard deviation close to 1.

Finally, if activation is not nothing, it is applied to the outputs as well.
"""
struct BatchNorm <: Layer

  	moments
	activation

	function BatchNorm(activation=identity)
		activation = get_activation(activation)
		new(Knet.bnmoments(), activation)
	end
end

function call(bn::BatchNorm, X::Tensor)
	bn.activation.(Knet.batchnorm(X, bn.moments))
end

@debug "Loaded Core modules"


text:
"""
### function plot_bode(gain, phase, w, num_u, num_y)

Plot the bode diagram.

##### Args

* gain: Gain [dB].
* phase: Phase [deg].
* w: Frequencies [rad/s].
* num_u: Number of inputs.
* num_y: Number of outputs.

"""

function plot_bode(gain, phase, w, num_u, num_y)
    if (plot_system == :pyplot)
        # Get the number of graphics.
        num_plots = length(gain[1])

        fig, ax = plt.subplots(num_y*2, num_u, sharex=true, sharey=true)

        for i = 1:num_plots
            gain_i  = map(x->x[i], gain)
            phase_i = map(x->x[i], phase)

            ax[2(i-1)+1][:semilogx](w, gain_i,  linewidth=plw)
            ax[2(i-1)+2][:semilogx](w, phase_i,  linewidth=plw)
            ax[2(i-1)+1][:grid]("on", which="major",
                         color="#AAAAAA", linewidth=glw)
            ax[2(i-1)+2][:grid]("on", which="major",
                         color="#AAAAAA", linewidth=glw)

            # Change figure border width.
            ax[2(i-1)+1][:spines]["top"][:set_linewidth](blw)
            ax[2(i-1)+1][:spines]["bottom"][:set_linewidth](blw)
            ax[2(i-1)+1][:spines]["left"][:set_linewidth](blw)
            ax[2(i-1)+1][:spines]["right"][:set_linewidth](blw)

            ax[2(i-1)+2][:spines]["top"][:set_linewidth](blw)
            ax[2(i-1)+2][:spines]["bottom"][:set_linewidth](blw)
            ax[2(i-1)+2][:spines]["left"][:set_linewidth](blw)
            ax[2(i-1)+2][:spines]["right"][:set_linewidth](blw)
        end

        for i = 1:num_u
            (num_u > 1) && ax[2*num_y*(i-1)+1][:set_title]("From Input $i")
            ax[2*num_y*i][:set_xlabel]("Freq. (rad/s)")
        end

        for i = 1:num_y
            ax[2(i-1)+1][:set_ylabel]("Gain (dB)")
            ax[2(i-1)+2][:set_ylabel]("Phase (deg)")

            ax_aux = ax[2(num_u-1)*num_y+2(i-1)+1][:twinx]()
            ax_aux[:yaxis][:set_label_position]("right")
            ax_aux[:set_yticks]([])
            ax_aux[:set_ylabel]("To Output $i")

            ax_aux = ax[2(num_u-1)*num_y+2(i-1)+2][:twinx]()
            ax_aux[:yaxis][:set_label_position]("right")
            ax_aux[:set_yticks]([])
            ax_aux[:set_ylabel]("To Output $i")
        end

        fig[:suptitle]("Bode Diagram", fontsize=tifsz)
    end
end


text:
function treshold_percent(array, percent=0.10; kwargs...)
    treshold =  floor(Int, length(array) * percent)
    return treshold == 0 ? 2 : treshold
end

function get_added_edges(clique1, clique2)
    edges = Set()
    for v1 in clique1
        for v2 in clique2
            if v1 != v2
                push!(edges, (v1, v2))
            end
        end
    end
    return edges
end

function merge_cliques(cliques, cliquetree, edge)
    c1 = cliques[edge[1]]
    c2 = cliques[edge[2]]
    new_cliques = filter(c -> !(c in [c1, c2]), deepcopy(cliques))
    push!(new_cliques, union(c1, c2))
    new_cliquetree = filter(e -> e != edge, deepcopy(cliquetree))
    update_cliquetree!(new_cliquetree, length(new_cliques), edge)
    return new_cliques, new_cliquetree
end

function update_cliquetree!(cliquetree, nb_cliques, edge)
    sorted_edge = sort(edge)
    for i in 1:length(cliquetree)
        for j in 1:2
            if cliquetree[i][j] in sorted_edge
                cliquetree[i][j] = nb_cliques
            elseif cliquetree[i][j] >= sorted_edge[2]
                cliquetree[i][j] -= 2
            elseif cliquetree[i][j] >= sorted_edge[1]
                cliquetree[i][j] -= 1
            end
        end
    end
end

# Used as a feature and not in the merge process
function compute_molzahn_heuristic(cliques, m, nb_lc)
    sizes = map(length, cliques)
    return m + nb_lc + sum(map(s -> s * ((s + 1) / 2), sizes))
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen