text:
# Type containing information on workers.
struct WorkerInf
    inputs::Dict{Int,RemoteChannel{Channel{Message}}}
    output::RemoteChannel{Channel{Message}}
    wait::Vector{Message}
    results::Dict{Int,Future}
end

function WorkerInf(inputs::Dict{Int,RemoteChannel{Channel{Message}}},
                   output::RemoteChannel{Channel{Message}},
                   results::Dict{Int,Future}
                   ) where T
    return WorkerInf(inputs,
                     output,
                     Vector{Message}(),
                     results)
end

function _abort(wi::WorkerInf)::Nothing
    for (wid, rc) in pairs(wi.inputs)
        put!(rc, Abort())
    end
    while isready(wi.output)
        # Ensure that workers aren't blocked
        take!(wi.output)
    end
    return
end

function _check_wait(wi::WorkerInf,
                     msg_type::Type{M}
                     )::Union{Nothing,M} where M <: Message

    # println("Checking wait queue...")

    retval = nothing

    for (k,msg) in enumerate(wi.wait)
        if typeof(msg) <: msg_type
            retval = msg
            deleteat!(wi.wait, k)
            break
        end
    end

    return retval
end

function _is_running(wi::WorkerInf)
    return length(wi.results) > 0
end

function _launch_workers(min_worker_size::Int,
                         min_result_size::Int,
                         )::WorkerInf

    worker_q_size = min_worker_size
    result_q_size = max(2*nworkers(), min_result_size)

    worker_results = RemoteChannel(()->Channel{Message}(result_q_size))
    worker_qs = Dict{Int,RemoteChannel{Channel{Message}}}()
    futures = Dict{Int,Future}()
    @sync for wid in workers()
        @async begin
            worker_q = RemoteChannel(()->Channel{Message}(worker_q_size), wid)
            futures[wid] = remotecall(worker_loop, wid, wid, worker_q, worker_results)
            worker_qs[wid] = worker_q
        end
    end

    wi = WorkerInf(worker_qs, worker_results, futures)
    _monitor_workers(wi)

    return wi
end

function _monitor_workers(wi::WorkerInf)::Nothing

    # println("Checking worker status...")

    for (pid, f) in pairs(wi.results)
        if isready(f)
            # Whatever the return value, this worker is done and
            # we can stop checking on it
            delete!(wi.results, pid)
            delete!(wi.inputs, pid)

            try

                # Fetch on a remote process result will throw and error here
                # whereas fetch on a task requires that it is thrown. Hence
                # the try-catch AND the type check for an exception.
                e = fetch(f)
                if typeof(e) <: Int
                    # Normal Exit -- do nothing
                elseif typeof(e) <: Exception
                    throw(e)
                else
                    error("Unexpected return code from worker: $e")
                end

            catch e

                _abort(wi)
                rethrow()

            finally

                ;

            end
        end
    end

    return

end

function _number_workers(wi::WorkerInf)
    return length(wi.results)
end

function _retrieve_message(wi::WorkerInf)::Message

    # println("Checking message queue...")

    while !isready(wi.output) && _is_running(wi)
        _monitor_workers(wi)
        yield()
    end

    msg = take!(wi.output)

    return msg
end

function _retrieve_message_type(wi::WorkerInf,
                                msg_type::Type{M}
                                )::M where M <: Message

    # Check if we have this type of message in the wait queue
    msg = _check_wait(wi, msg_type)

    while isnothing(msg)

        # Check get messages from the remote channel
        msg = _retrieve_message(wi)

        # If it's wrong put it on the wait queue for later and keep looking
        if !(typeof(msg) <: msg_type)
            push!(wi.wait, msg)
            msg = nothing
        end

    end

    return msg

end

function _send_message(wi::WorkerInf,
                       wid::Int,
                       msg::Message,
                       )
    # println("Sending message of type $(typeof(msg)) to worker $(wid)")
    put!(wi.inputs[wid], msg)
    return
end

function _shutdown(wi::WorkerInf)
    @sync for rchan in values(wi.inputs)
        @async put!(rchan, ShutDown())
    end
    return
end

function _wait_for_shutdown(wi::WorkerInf)::Nothing
    while _is_running(wi)
        _monitor_workers(wi)
        yield()
    end
    return
end


text:
module sampler

import StatsBase
using ergm.spaces, ergm.models
import ergm.stats

export GibbsSampler, sample, ParallelGibbsSampler

mutable struct GibbsSampler
    initial_state
    state
    model
    burn_in
    sample_interval

    function GibbsSampler(
            initial_state,
            model :: ExponentialFamily,
            burn_in,
            sample_interval
    )
        sampler = new(initial_state, undef, model, burn_in, sample_interval)
        restart(sampler)
        sampler
    end
end

function restart(sampler :: GibbsSampler)
    sampler.state = copy(sampler.initial_state)
    set_state(sampler.model, sampler.state)
end

function gibbs_step(sampler :: GibbsSampler)
    is = keys(sampler.state)
    n = length(is)

    for i ∈ is
        d = getdomain(sampler.state, i)
        w = zeros(length(d))
        old_x = sampler.state[i]

        for (j, x) ∈ enumerate(d)
            w[j] = test_state(sampler.model, (i, x))
        end

        # Try to avoid going outside Float64
        # precision when applying exp. Can
        # exponentiation be avoided here? Can
        # be avoided for length(d) == 2 but not
        # sure how to generalize.
        c = (maximum(w) + minimum(w)) / 2
        w = exp.(w .- c)
        x = StatsBase.sample(d, StatsBase.Weights(w))
        sampler.state[i] = x
        update_state(sampler.model, (i, x))
    end
end

function sample(sampler :: GibbsSampler, n)
    restart(sampler)
    samples = []
    m = length(sampler.model.params)
    sample_stats = zeros(m, n)

    # burn in to reach equilibrium state
    for _ ∈ 1:sampler.burn_in
        gibbs_step(sampler)
    end

    for i ∈ 1:n
        # draw one sample
        push!(samples, copy(sampler.state))
        sample_stats[:, i] = stats.get_stats(sampler.model.stats)

        # throw away some samples to reduce autocorrelation
        for _ ∈ sampler.sample_interval
            gibbs_step(sampler)
        end
    end

    samples, sample_stats
end

struct ParallelGibbsSampler
    n_samplers
    samplers :: Vector{GibbsSampler}

    function ParallelGibbsSampler(
            initial_state,
            model :: ExponentialFamily,
            burn_in,
            sample_interval,
            n_samplers
    )
        samplers = [
            GibbsSampler(initial_state, copy(model), burn_in, sample_interval)
            for _ ∈ 1:n_samplers
        ]
        new(n_samplers, samplers)
    end
end

function sample(sampler :: ParallelGibbsSampler, n)
    # split requested samples evenly across chains
    m, r = divrem(n, sampler.n_samplers)
    n_per_chain = vcat(
        fill(m + 1, r),
        fill(m, sampler.n_samplers - r)
    )
    hs = []

    for i ∈ 1:sampler.n_samplers
        h = Threads.@spawn sample(sampler.samplers[i], n_per_chain[i])
        push!(hs, h)
    end

    np = length(get_params(sampler.samplers[1].model))
    samples = []
    sample_stats = zeros(np, n)
    j = 1

    for (i, h) ∈ enumerate(hs)
        s, ss = fetch(h)
        samples = vcat(samples, s)
        sample_stats[:, j:j + n_per_chain[i] - 1] = ss
        j += n_per_chain[i]
    end

    samples, sample_stats
end

end

ext: jl


text:
using Base: @propagate_inbounds

eltype_extension(a::AbstractArray) = eltype(a)

@inline index_extension(_, _, ext::Extension) =
    throw(ArgumentError("index_extension() undefined for extension $ext"))

@inline index_extension(a, i, ext::ExtensionReplicate) =
    CartesianIndex(clamp.(Tuple(i), 1, size(a)))

@inline index_extension(a, i, ext::ExtensionCircular) =
    CartesianIndex(mod1.(Tuple(i), size(a)))

symidx(k, n) = mod1(isodd(fld1(k, n)) ? k : -k, n)
@inline index_extension(a, i, ext::ExtensionSymmetric) =
    CartesianIndex(symidx.(Tuple(i), size(a)))

"""
    getindex_extension(w::Window, wi::CartesianCoordinate, ext::Extension)

Return a value for an indexing operation on `w` with relative index `wi`
respecting the extension behaviour `ext`.
Users may define their own extension behaviour by defining additional methods
for this function.
"""
getindex_extension

"""
    setindex_extension!(w::Window, x, wi::CartesianCoordinate, ext::Extension)

Store the value `x` for an indexing operation on `w` with relative index `wi`
respecting the extension behaviour `ext`.
Users may define their own extension behaviour by defining additional methods
for this function.
Depending on the extension a call to this function may be non-sensical and thus
may be implemented as a no-op.
"""
setindex_extension!

@propagate_inbounds getindex_extension(w::Window, wi) =
    getindex_extension(parent(w), position(w) + wi)

@propagate_inbounds setindex_extension!(w::Window, x, wi) =
    setindex_extension!(parent(w), x, position(w) + wi)

@inline getindex_extension(_, _, ext::ExtensionNothing) = nothing
@inline setindex_extension!(_, _, _, ext::ExtensionNothing) = nothing

@inline getindex_extension(_, _, ext::ExtensionConstant) = ext.value
@inline setindex_extension!(_, _, _, ext::ExtensionConstant) = ext.value

@propagate_inbounds getindex_extension(a, i, _) =
    a[index_extension(a, i, extension(a))]

@propagate_inbounds setindex_extension!(a, x, i, _) =
    a[index_extension(a, i, extension(a))] = x


text:
##############################################################################
#
# LDA
#
##############################################################################

lda(dtm::DocumentTermMatrix) = error("not yet implemented")
lda(crps::Corpus) = error("not yet implemented")


text:
using DICOM, Lucas
using DocOpt

args = docopt("""showdicom.jl

Usage:
  dicom.jl show [-v] [--tag=<tag>] (--dir=<dir> | <files>...)
  dicom.jl (dump | debug) [-v] (--dir=<dir> | <files>...)
  dicom.jl -h | --help

Options:
  -h --help     Show this screen.
  -v --verbose  Show more verbose output.
  --tag=<tag>   Shows only one specific tag (tag syntax is 0xabcd,0x1234).
  --dir=<dir>   Use all files of `dir`, which is either `train` or `test`.
""")

if args["--tag"] != nothing
    args["--tag"] = tuple([parse(UInt16, t) for t in split(args["--tag"], ',')]...)
end

if args["--dir"] != nothing
    println("Collecting files...")
    args["<files>"] = UTF8String[]
    for d in readdir(args["--dir"])
        patientdir = joinpath(args["--dir"], d, "study")
        for s in readdir(patientdir)
            for f in readdir(joinpath(patientdir, s))
                if splitext(f)[2] == ".dcm"
                    push!(args["<files>"], joinpath(patientdir, s, f))
                end
            end
        end
    end
end

for name in args["<files>"]
    if args["--verbose"]
        println("FILE: ", name)
        println("=====")
    end
    dcm = dcm_parse(open(name));

    if args["show"]
        for el in dcm
            if args["--tag"] == nothing || el.tag == args["--tag"]
                print(el.tag, ": ")
                display(el)
                println()
            end
        end
    elseif args["dump"]
        # Dumps the image into a csv file.
        if (el = DICOM.lookup(dcm, (0x7fe0,0x0010))) != false
            writecsv("$name.csv", el.data[1][:,:,1])
            # TODO: Possibly also dump these. at least look into them!
            tags = Dict(
                (0x0008,0x0008) => "imtype",  # Image Type
                (0x0008,0x0060) => "modality",  # Modality: ASCIIString["MR"]
                (0x0008,0x0070) => "manuf",  # Manufacturer: ASCIIString["SIEMENS"]
                (0x0008,0x103e) => "series",  # Series Description: ASCIIString["2ch"]
                (0x0008,0x1090) => "model",  # Manufacturer's Model Name: ASCIIString["Aera"]
                (0x0018,0x0022) => "opts",  # Scan Options: ASCIIString["CT"]
                (0x0018,0x0023) => "aqtype",  # MR Acquisition Type: ASCIIString["2D"]
                (0x0018,0x0050) => "slice_thick",  # Slice Thickness: [8.0]
                (0x0018,0x1314) => "flip_angle",  # Flip Angle: [56.0]
                (0x0018,0x5100) => "pos",  # Patient Position: ASCIIString["HFS"]
                (0x0020,0x0032) => "imgpos",  # Image Position (Patient): [-87.771005468217,93.799321675558,195.08339402037]
                (0x0020,0x0037) => "imgori",  # Image Orientation (Patient): [0.8127451738409,-0.5549969637339,0.1772671787004,0.2130993502411,3.956532e-8,-0.9770305353093]
                (0x0020,0x1041) => "sliceloc",  # Slice Location: [53.505870937259]
                (0x0028,0x0004) => "photometric",  # Photometric Interpretation: ASCIIString["MONOCHROME2"]
                (0x0028,0x0030) => "pixelsp",  # Pixel Spacing: [1.40625,1.40625]
                (0x0028,0x1050) => "win_cent",  # Window Center: [344.0]
                (0x0028,0x1051) => "win_width",  # Window Width: [848.0]
                (0x0028,0x1055) => "win_expl",  # Window Center & Width Explanation: ASCIIString["Algo1"]
                (0x7fe0,0x0010) => "data",  # Variable Pixel Data: 256x192x1 Array{UInt16,3}:
            )
        end
    elseif args["debug"]
        # TODO: Write whatever code you're curious about right now here!

        # Show image size
        #if (el = DICOM.lookup(dcm, (0x7fe0,0x0010))) != false
        #    println(size(el.data[1]))
        #end

        # Check if "max pixel" == max(pixels) or not
        # Yes, it is.
        #maxpix = DICOM.lookup(dcm, (0x0028,0x0107)).data[1]
        #imgdata = DICOM.lookup(dcm, (0x7fe0,0x0010)).data[1]
        #println(maximum(imgdata) == maxpix)
    end
end


text:
# TODO keep track of outgoing arcs?
"""
    DependencyTree

A rooted tree of dependency relations among the tokens of a sentence.
"""
struct DependencyTree{T<:Token, R<:Union{Int,Set{Int}}}
    tokens::Vector{T}
    root::R
end

"""
    deptree(read_token, xs)

Create a `DependencyTree` by calling `read_token` on each of `xs`.
"""
function deptree end

function deptree(tokens)
    return DependencyTree(tokens, find_root(tokens))
end

function deptree(read_token, xs)
    tokens = Token[]
    for x in xs
        try
            token = read_token(x)
            push!(tokens, token)
        catch err
            if err isa EmptyTokenError
                continue
            elseif err isa MultiWordTokenError
                continue
            else
                throw(err)
            end
        end
    end
    tokens = identity.(tokens) # type info
    root = find_root(tokens)
    return DependencyTree(tokens, root)
end

function deptree(lines::String, read_token=from_conllu)
    lines = [line for line in split(strip(lines), "\n")]
    lines = filter(x -> !startswith(x, "#"), lines)
    return deptree(read_token, lines)
end

function find_root(tokens)
    root = [i for (i, token) in enumerate(tokens) if has_head(token, 0)]
    return length(root) == 1 ? root[1] : Set(root)
end

"""
    arcs(tree::DependencyTree)

Return a vector of (src, dst) tuples, representing all dependency arcs in `tree`.
"""
arcs(tree::DependencyTree) =
    [(h, d) for (d, tok) in enumerate(tree.tokens) for h in tok.head]

deps(tree::DependencyTree, i::Int) =
    filter(j -> has_head(tree.tokens[j], i), 1:length(tree))

leftdeps(tree::DependencyTree, i::Int) = filter(j -> j < i, deps(tree, i))

function leftmostdep(tree::DependencyTree, i::Int)
    ds = leftdeps(tree, i)
    isempty(ds) ? -1 : first(ds)
end

rightdeps(tree::DependencyTree, i::Int) = filter(j -> j > i, deps(tree, i))

function rightmostdep(tree::DependencyTree, i::Int)
    ds = rightdeps(tree, i)
    isempty(ds) ? -1 : last(ds)
end

has_arc(tree::DependencyTree, head::Int, dependant::Int) =
    iszero(dependant) ? false : has_head(tree.tokens[dependant], head)

children(tree::DependencyTree; recursive=false) =
    [children(tree, i; recursive=recursive) for i in 1:length(tree)]

function children(tree::DependencyTree, i::Int; recursive=false)
    ch = [c for (c, tok) in enumerate(tree.tokens) if has_head(tok, i)]
    if recursive
        for c in ch
            append!(ch, children(tree, c))
        end
    end
    return ch
end

token(tree::DependencyTree, i) =
    iszero(i) ? ROOT : tree.tokens[i]

function is_projective(tree::DependencyTree)
    # For every arc (i,l,j) there is a directed path from i to every
    # word k such that min(i,j) < k < max(i,j)
    arc = (a, b) -> has_arc(tree, a, b)
    for (parent, child) in arcs(tree)
        child > parent && ((child, parent) = (parent, child))
        for k in child+1:parent-1
            for m in 1:length(tree)#Iterators.flatten((1:i-1, j+1:length(tree)))
                if m < child || m > parent
                    if arc(k, m) || arc(m, k)
                        return false
                    end
                end
            end
        end
    end
    return true
end

struct ArrowCharset
    arrow::Char
    ur::Char # up/right
    ud::Char # up/down
    lr::Char # left/right
    dr::Char # down/right
end
DEFAULT_ARROWS = ArrowCharset('►','└','│','─','┌')
ASCII_ARROWS   = ArrowCharset('>', '\'', '|', '-', ',')

function prettyprint(tree::DependencyTree; charset=DEFAULT_ARROWS)
    n = length(tree)
    lines = ["" for _=1:n+1]
    queue = sort(arcs(tree), by = (arc) -> abs(-(arc...)))
    for (h, d) in queue
        iszero(h) && continue # leave root arc until the end
        lines = drawarrow!(lines, h+1, d+1, charset)
    end
    for root in tree.root
        lines = drawarrow!(lines, 1, root+1, charset)
    end
    max_height = maximum(length.(lines))
    lines = map(enumerate(lines)) do (i, line)
        height = length(line)
        leftpad = repeat(" ", max_height - height)
        leftpad * line * " " * tree[i-1].form
    end
    return join(lines, "\n")
end

function drawarrow!(lines, from, to, cs)
    middle = cs.ud * "  "
    if from < to # pointing "down" the sentence
        range = from:to
        base = cs.dr * cs.lr * cs.lr
        head = cs.ur * cs.lr * cs.arrow
        arrow = [base; [middle for i in from+1:to-1]; head]
    else
        range = to:from
        arrow = [cs.dr * cs.lr * cs.arrow;
                 [middle for i in to+1:from-1];
                 cs.ur * cs.lr * cs.lr]
    end
    heights = length.(lines[i] for i in range)
    max_height = maximum(heights)
    for (i, (r, h, a)) in enumerate(zip(range, heights, arrow))
        if r == from
            lines[r] = a * repeat(cs.lr, max_height - h) * lines[r]
        elseif r == to
            # redraw the arrowhead to go all the way down
            arr = a[1][1] * repeat(cs.lr, max(0, max_height - h) + 1) * cs.arrow
            lines[r] = arr * lines[r]
        else
            lines[r] = a * repeat(" ", max_height - h) * lines[r]
        end
    end
    return lines
end

Base.getindex(tree::DependencyTree, i::Int) = iszero(i) ? ROOT : tree.tokens[i]

Base.iterate(tree::DependencyTree, state...) = iterate(tree.tokens, state...)

Base.length(tree::DependencyTree) = length(tree.tokens)

Base.show(io::IO, tree::DependencyTree) =
    print(io, prettyprint(tree))

==(a::DependencyTree, b::DependencyTree) = all(a.tokens .== b.tokens)


text:
module NestedPseudoLikelihood
using StructuralBase, MDPTools, StatsFuns, Optim
import StructuralBase: inner_iterations, fit
export fit_npl, fit

include("types.jl")
include("logit.jl")
include("initial_ccps.jl")

# The StatsBase conformant way of writing it
fit(E::NPL, U, X, D) = fit_npl(U, X, D, E)

# In case someone passes a single state
fit_npl(U, X::MDPTools.DiscreteState, D; kwargs...) = fit_npl(U, States(X), D; kwargs...)

# The
function fit_npl(U, X, D; method = Optim.Newton())
    E = NPL(U, fit_np!(U.P, D.a, D.x))
    fit_npl(U, X, D, E; method = method)
end
function fit_npl{Ts<:MDPTools.AbstractState}(U::LinearUtility,
                                            X::Ts,
                                            D::Data,
                                            E::NPL;
                                            ## Keyword arguments
                                            # Convergence parameters
                                            ε_outer = 1e-8, # ununsed
                                            ε_inner = 1e-8, # unused
                                            # Initial estimation
                                            method::Optim.Optimizer = Optim.Newton())
    # Get dimension of choice vector once and for all
    J = length(U.P)

    # Initialize cache variables such that most calculations are in-place
    c = NPLCacheVars(X.nX, length(D.xᵈ), U.nvar, J, X.F)

    # Keep track of convergence info
    conv_info = ConvergenceNPL()

    # Keep a trace of parameters and log-likelihood values
    tr = TraceNPL(E.K, U.nvar)

    # Define short-hand notation for general likelihood expressions at data
        llᵈ(x)          =     ll(x,          D, c, J)
       ∇llᵈ(x, grad)    =    ∇ll(x, grad,    D, c, J)
    ll_∇llᵈ(x, grad)    = ll_∇ll(x, grad,    D, c, J)
      ∇²llᵈ(x, hessian) =   ∇²ll(x, hessian, D, c, J)

    # Set options; could be a positional argument?
    options = Optim.Options(show_trace = E.verbose, extended_trace=E.verbose)

    for k = 1:E.K
        copy!(c.Θᵏ, U.Θ)

        ## Prepare inner nest
        # Solve for WP given current P
        Wᴾ!(c, U, X, J)

        # Update the variables used to construct the choice specific values
        discounted_vars!(c, U, X, J)
        # Update variables used in the conditional likelihood
          observed_vars!(D, c, J)

        # Optimize likelihood and apply policy update
        likelihood = TwiceDifferentiable(llᵈ, ∇llᵈ, ll_∇llᵈ, ∇²llᵈ, c.Θᵏ) # sort of a hack, could be waste
        res = optimize(likelihood, c.Θᵏ, method, options)

        # trace the results from the maximum likelihood step
        trace!(tr, conv_info, res, U, E, k)

        # Perform policy step given ᵏ variables to get Pᵏ⁺¹
        Ψ!(c, U)

        # If we're not estimating NPL-K, assess convergence
        !E.finite && assess_convergence(norm(c.Θᵏ - U.Θ, Inf), ε_outer, E, conv_info, k)
    end

    # Store final gradient and Hessian for users to inspect
    gradient = zeros(U.nvar)
    hessian = zeros(U.nvar, U.nvar)
    _ll = ll_∇llᵈ(c.Θᵏ, gradient)
            ∇²llᵈ(c.Θᵏ, hessian)

    # Return results
    objective = TwiceDifferentiable(llᵈ, ∇llᵈ, ll_∇llᵈ, ∇²llᵈ, c.Θᵏ)
    EstimationResults(E,
                      -_ll*D.nobs,
                      -gradient*D.nobs,
                      -hessian*D.nobs,
                      objective,
                      copy(U.Θ),
                      conv_info,
                      tr,
                      D.nobs,
                      objective)
end

"""
    Wᴾ!(c, U, X, J)
    Wᴾ!(  Fᵁ,   Wᴾ,   β,   F,   eᴾ,   Pᵏ,   Z, J)
Updates the matrix Wᴾ that contains the variables used to calculate the valuation
given a current Pᵏ and parameters Θᵏ, Wᴾz*Θᵏ+Wᵖe.
"""
Wᴾ!(c, U, X, J) =
         Wᴾ!(c.Fᵁ, c.Wᴾ, U.β, X.F, c.eᴾ, U.P, U.Z, J)
function Wᴾ!(  Fᵁ,   Wᴾ,   β,   F,   eᴾ,   Pᵏ,   Z, J)
    # Calculate E(ϵʲ|j optimal)
    map!(p->-log.(p), eᴾ, Pᵏ) # this should dispatch on some shock type

    # Calculate transition probability matrix given current policy Pᵏ
    Fᵁ .= β*sum(Pᵏ[ia].*F[ia] for ia = 1:J)

    # Update Wᴾ = [Wᵖz Wᵖe]
    Wᴾ .= (I-Fᵁ)\sum(Pᵏ[ia].*[Z[ia] eᴾ[ia]] for ia = 1:J)
end

discounted_vars!(c, U, X, J) =
         discounted_vars!(c.z̃, c.ẽ, c.Δz̃, c.Δẽ, c.Wᴾ, U, X, J)
function discounted_vars!(  z̃,   ẽ,   Δz̃,   Δẽ,   Wᴾ, U, X, J)
    @views for j = J:-1:1
        z̃[j] .= U.Z[j] + U.β * (X.F[j] * Wᴾ[:,1:end-1])
        ẽ[j] .= U.β * X.F[j] * Wᴾ[:,end]
        if j < J
            Δz̃[j] .= z̃[j]-z̃[J]
            Δẽ[j] .= ẽ[j]-ẽ[J]
        end
    end
end

function observed_vars!(D, c, J)
    # Only use variables at the observed states
    for j = 1:J
      c.z̃ᵈ[j] .= c.z̃[j][D.xᵈ,:]
        if j < J
            c.Δz̃ᵈ[j] .= c.Δz̃[j][D.xᵈ,:]
            c.Δẽᵈ[j] .= c.Δẽ[j][D.xᵈ]
        end
    end
end

"""
    Ψ!(c, U)
    Ψ!(Δz̃, Δẽ, m, U, P)

The best response operator that takes ᵏ variables and updates Pᵏ to Pᵏ⁺¹ in-place.
"""
Ψ!(c, U) =
         Ψ!(c.Δz̃, c.Δẽ, c.m, U, U.P)
# Note that we're updating the whole state space, not just xᵈ, as we need Pᵏ⁺¹
# at all x ∈ X to construct Fᵁ that is needed to calculate Wᴾ.
function Ψ!(Δz̃, Δẽ, m, U, P)
    J = length(P)
    for i = 1:J-1
        m[i] .= Δz̃[i]*U.Θ+Δẽ[i]
    end
    denom = 1+sum(exp.(m[i]) for i = 1:J-1)
    for i = 1:J-1
        P[i] .= exp.(m[i])./denom
    end
    P[J] .= 1-sum(P[ia] for ia = 1:J-1)
end

function assess_convergence(norm, ε_outer, E, conv_info, k)
    if norm < ε_outer
        conv_info.flag = true
        conv_info.outer = k
        conv_info.norm_Θ = norm
        E.Θu = E.Θu[1:k, :]
        E.ll = E.ll[1:k]
    end
end

function trace!(tr, conv_info, res, U, E, k)
    # Increment iteration counter
    conv_info.iter_maxlike += Optim.iterations(res)

    # Grab estimates
    copy!(U.Θ, Optim.minimizer(res))

    # Trace loglikelihood
    E.ll[k] = conv_info.ll

    # Trace parameter estimates
    E.Θu[k, :] = U.Θ
    tr.Θ[k,:] = U.Θ
end

end # module


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, "..")); Pkg.instantiate()
Pkg.activate(); Pkg.instantiate()

pushfirst!(LOAD_PATH, joinpath(@__DIR__, ".."))

using Documenter, SourceWalk

makedocs(
    sitename = "SourceWalk",
    pages = [
        "Home" => "index.md"],
    format = Documenter.HTML(prettyurls = haskey(ENV, "CI")))

deploydocs(
  repo = "github.com/MikeInnes/SourceWalk.jl.git",)


text:
function  neightriplets(x, q, tripletlist, tripletnum, atomtype, ilist, alist)

    dim = size(x,1);
    ncq = size(q,1);
    inum = length(ilist);
    ijknum = sum(tripletnum[ilist.+1]-tripletnum[ilist]);

    ai = zeros(Int64,ijknum);
    aj = zeros(Int64,ijknum);
    ak = zeros(Int64,ijknum);
    ti = zeros(Int64,ijknum);
    tj = zeros(Int64,ijknum);
    tk = zeros(Int64,ijknum);
    xij = zeros(dim,ijknum);
    xik = zeros(dim,ijknum);
    qi = zeros(ncq,ijknum);
    qj = zeros(ncq,ijknum);
    qk = zeros(ncq,ijknum);

    count = 0;
    for ii=1:inum
        i = ilist[ii]; # atom i
        n1 = tripletnum[i];
        n2 = tripletnum[i+1];
        n = n2-n1;
        gj = tripletlist[(n1+1):n2,1]; # a list of neighbors around atom i
        j = alist[gj];
        gk = tripletlist[(n1+1):n2,2]; # a list of neighbors around atom i
        k = alist[gk];

        ind = (count+1):(count+n);
        count = count + n;

        ai[ind] .= i;
        aj[ind] = j;
        ak[ind] = k;
        ti[ind] .= atomtype[i];
        tj[ind] = atomtype[j];
        tk[ind] = atomtype[k];
        xij[:,ind] = x[:,gj] .- x[:,i];
        xik[:,ind] = x[:,gk] .- x[:,i];

        # display(i)
        # display(j)
        # display(k)
        # display(ind)
        # display(xij[:,ind])
        # display(xik[:,ind])
        # rij = sqrt.(xij[1,ind].*xij[1,ind] + xij[2,ind].*xij[2,ind] + xij[3,ind].*xij[3,ind])
        # display(maximum(rij))
        # rij = sqrt.(xik[1,ind].*xik[1,ind] + xik[2,ind].*xik[2,ind] + xik[3,ind].*xik[3,ind])
        # display(maximum(rij))
        # error("here")
        if ncq>0
            qi[ind] .= q[i];
            qj[ind] = q[j];
            qk[ind] = q[k];
        end
    end

    return xij, xik, ai, aj, ak, ti, tj, tk, qi, qj, qk

    end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
# # Two Parabolas
#
# The “two parabolas” problem in two dimensions reads as
# ```math
#     \min_{𝐱 ∈ X }
#     \begin{bmatrix} f₁(\mathbf{x}) \\ f₂(\mathbf{x}) \end{bmatrix} =
#     \min_{\mathbf{x} ∈ X}
#     \begin{bmatrix}
#     (x₁ - 1)² + (x₂ - 1)² \\
#     (x₁ + 1)² + (x₂ + 1)²
#     \end{bmatrix}.
# ```
# It is unconstrained if the feasible set is ``X = ℝ^2``.
# The individual minima ``[1,1]`` and ``[-1,-1]`` are such that (in the unconstrained case)
# the global Pareto Set is
# ```math
# \mathcal{P}_{S} = \{ \mathbf{x} ∈ ℝ^2 : x₁ = x₂, \, -1 \le x₁, x₂ \le 1  \}.
# ```
#
# ## Solve using Exact Functions
#
# The gradients are easily calculated as
# ```math
# \nabla f_1 (\mathbf x) = 2 \begin{bmatrix}
# x_1 -1 \\ x_2 - 1 \end{bmatrix}, \;
# \nabla f_2 (\mathbf x) = 2 \begin{bmatrix}
# x_1 +1 \\ x_2 + 1 \end{bmatrix}, \;
# ```
#
# We can provide them to the solver to find a critical point:

using Pkg #src
Pkg.activate(@__DIR__) #src
using Test #src

using Morbit
Morbit.print_all_logs() #src

f₁ = x -> sum( (x .- 1).^2 )
f₂ = x -> sum( (x .+ 1).^2 )
∇f₁ = x -> 2 .* ( x .- 1 )
∇f₂ = x -> 2 .* ( x .+ 1 )

mop = MixedMOP(2);  # problem with 2 variables
add_objective!(mop, f₁, ∇f₁ )
add_objective!(mop, f₂, ∇f₂ )

#~ starting point
x₀ = [ -π ;  2.71828 ]

#~ set maximum number of iterations
ac = AlgoConfig( max_iter = 20)
#~ `optimize` will return parameter and result vectors as well
#~ as an return code and the evaluation database:
x, fx, ret_code, db = optimize( mop, x₀; algo_config = ac );
x

# Hopefully, `x` is critical, i.e., `x[1] ≈ x[2]`.
@test x[1] ≈ x[2] atol = .1 #src

# !!! note
#     To print more information on what the solver is doing, you can use the `Logging` module:
#     ```julia
#     import Logging: global_logger, ConsoleLogger
#     global_logger( ConsoleLogger( stderr, Morbit.loglevel4;
#         meta_formatter = Morbit.morbit_formatter ) )
#     ```
#     `loglevel4` is the most detailed and `loglevel1` is least detailed.
#     `Morbit.print_all_logs()` is a convenient shorthand.

#%% #src
# ### Plotting Iteration Sites
# Let's retrieve the iteration sites.
# We convert to Tuples for easier plotting.
iteration_indices = [ iter_.x_index for iter_ in db.iter_info]
it_sites = Tuple.(Morbit.get_site.(db, iteration_indices))

# For Plotting we use CairoMakie
using Makie, CairoMakie

#~ Pareto Set ≙ line from (-1,-1) to (1,1)
fig, ax, _ = lines( [(-1,-1),(1,1)]; color = :blue, linewidth = 2,
    figure = (resolution = (600, 600),) )

#~ Plot the iteration sites:
lines!(it_sites)
scatter!(it_sites;
    color = LinRange(0, 1, length(it_sites)),
    colormap = :winter
)

#~ Plot function contours
Y = X = LinRange(-4, 4, 100)
Z₁ = [ f₁([x;y]) for x ∈ X, y ∈ X ]
Z₂ = [ f₂([x;y]) for x ∈ X, y ∈ X ]
levels = [ i.^2 for i = LinRange(.1, 6, 6) ]
contour!(X,Y,Z₁; colormap = :greens, levels = levels, linewidth = .5 )
contour!(X,Y,Z₂; colormap = :heat, levels = levels, linewidth = .5 )

#~ Show the plot:
ax.title[] = "Pareto Set and Iterates."
ax.xgridvisible[] = false
ax.ygridvisible[] = false

fig

#%% #src

# ## Solving using RBF Surrogates
#
# Suppose now that we do not have access to the objective gradients and that the objectives
# also take some time to evaluate.
# In this situation, we could try to model them using surrogate models.
# To use radial basis function models, pass an `RbfConfig` when specifying the objective:
#
mop_rbf = MixedMOP()

#~ Define the RBF surrogates
rbf_cfg = RbfConfig(
    kernel = :inv_multiquadric
)
#~ Add objective functions to `mop_rbf`
add_objective!(mop_rbf, f₁, rbf_cfg )
add_objective!(mop_rbf, f₂, rbf_cfg )

#~ only perform 10 iterations
ac = AlgoConfig( max_iter = 10 )
x, fx, _, db = optimize( mop_rbf, x₀; algo_config = ac )
x

#src Setup and save plot for docs.
iteration_indices_rbf = [ iter_.x_index for iter_ in db.iter_info]
it_sites_rbf = Tuple.(Morbit.get_site.(db, iteration_indices_rbf))
lines!(it_sites_rbf) #hide
scatter!(it_sites_rbf; color = :orange) #hide
nothing #hide

# The iteration sites are the orange circles:
fig #hide

#%% #src

# ## Different Starting Points and Recycling Data
#
# The method could converge to different points depending on the starting point.
# We can pass the evaluation data from previous runs to facilitate the construction of surrogate models:

#src Setup the problem anew, to ensure fresh start
ac = AlgoConfig( #hide
    max_iter = 10 #hide
    ); #hide
mop_rbf = MixedMOP(); #hide
#~ define the RBF surogates #hide
rbf_cfg = RbfConfig(  #hide
    kernel = :inv_multiquadric, #hide
); #hide
#~ add objective functions to `mop_rbf` #hide
add_objective!(mop_rbf, f₁, rbf_cfg ); #hide
add_objective!(mop_rbf, f₂, rbf_cfg ); #hide

#~ an array of well spread points in [-4,4]² #hide
X =[ #hide
 [-4.0, -4.0], #hide
 [3.727327839472812, 3.8615291196035457], #hide
 [3.804712690019901, -3.9610212058521235], #hide
 [-0.14512898384374573, -0.005775390168885508], #hide
 [-3.775315499879552, 3.8150054323309064], #hide
 [1.714228746087743, 1.8435786475209621], #hide
 [-1.9603720505875337, -2.0123206708499275], #hide
 [3.9953803225349187, -0.47734576293976794], #hide
 [-3.9944468955728745, 0.49857343385493635], #hide
 [-1.0455585089057458, 2.735699160002545] #hide
]; #hide

# Suppose, `X` is a list of different points in ℝ².

#src This is the code block visible in the docs:

#~ A dict to associate starting and end points:
start_fin_points = Dict();

#~ perform several runs:
db₀ = nothing # initial database can be `nothing`
for x₀ ∈ X
    global db₀, start_fin_points
    x_fin, fx_fin, _, db₀ = optimize( mop_rbf, x₀; algo_config = ac, populated_db = db₀ )
    #~ add points to dict
    start_fin_points[x₀] = x_fin
end

# Plotting:

fig, ax, _ = lines( [(-1,-1),(1,1)]; color = :blue, linewidth = 2,
    figure = (resolution = (600, 600), ),
    axis = (title="Different Starting Points",),
)

for (k,v) in start_fin_points
    lines!( [ Tuple(k), Tuple(v) ]; color = :lightgray )
end

scatter!( Tuple.(keys(start_fin_points));
    color = :green
)
scatter!( Tuple.(values(start_fin_points));
    color = :lightblue
)

fig #hide

# In the plot, the green points show the starting points and the lightblue circles show the final iterates:


text:
abstract Optimizer

# NLoptimizer holds the parameters passed to Opt. This is a workaround for the lack of serialization of Opt.
immutable NLoptimizer <: Optimizer
  npars::Int # Number of optimization parameters
  alg::Symbol # Optimization algorithm to use
  lb::VF64OrNothing # Lower bound constraint
  ub::VF64OrNothing # Upper bound constraint
  absftol::Float64 # Stopping criterion: absolute tolerance on function value
  relftol::Float64 # Stopping criterion: relative tolerance on function value
  absxtol::Float64 # Stopping criterion: absolute tolerance on optimization parameters
  relxtol::Float64 # Stopping criterion: relative tolerance on optimization parameters
  nfeval::Int64 # Stopping criterion: maximum number of function evaluations
  maxtime::Float64 # Stopping criterion: maximum optimization time

  function NLoptimizer(npars::Int, alg::Symbol, lb::VF64OrNothing, ub::VF64OrNothing, absftol::Float64,
    relftol::Float64, absxtol::Float64, relxtol::Float64, nfeval::Int64, maxtime::Float64)
    if isnan(absftol) && isnan(relftol) && isnan(absxtol) && isnan(relxtol) && nfeval < 0 && maxtime < 0
      error("A stopping criterion must be specified.")
    end

    new(npars, alg, lb, ub, absftol, relftol, absxtol, relxtol, nfeval, maxtime)
  end
end

NLoptimizer(npars::Int;
  alg::Symbol=:LN_COBYLA,
  lb::VF64OrNothing=nothing,
  ub::VF64OrNothing=nothing,
  absftol::Float64=1e-32,
  relftol::Float64=NaN,
  absxtol::Float64=1e-32,
  relxtol::Float64=NaN,
  nfeval::Int64=1_000,
  maxtime::Float64=-1.0) =
  NLoptimizer(npars, alg, lb, ub, absftol, relftol, absxtol, relxtol, nfeval, maxtime)

function NLoptimizer(;
  alg::Symbol=:LN_COBYLA,
  lb::F64OrNothing=nothing,
  ub::F64OrNothing=nothing,
  absftol::Float64=1e-32,
  relftol::Float64=NaN,
  absxtol::Float64=1e-32,
  relxtol::Float64=NaN,
  nfeval::Int64=1_000,
  maxtime::Float64=-1.0)
  lbvec = lb == nothing ? nothing : [lb]
  ubvec = ub == nothing ? nothing : [ub]
  NLoptimizer(1, alg, lbvec, ubvec, absftol, relftol, absxtol, relxtol, nfeval, maxtime)
end

function convert(::Type{Opt}, o::NLoptimizer)
  opt = Opt(o.alg, o.npars)

  if o.lb != nothing; lower_bounds!(opt, o.lb); end
  if o.ub != nothing; upper_bounds!(opt, o.ub); end
  if !isnan(o.absftol) ftol_abs!(opt, o.absftol); end
  if !isnan(o.relftol) ftol_rel!(opt, o.relftol); end
  if !isnan(o.absxtol) xtol_abs!(opt, o.absxtol); end
  if !isnan(o.relxtol) xtol_rel!(opt, o.relxtol); end
  if !(o.nfeval < 0) maxeval!(opt, o.nfeval); end
  if !(o.maxtime < 0) maxtime!(opt, o.maxtime); end

  return opt
end


text:
module Latent

include("GMM.jl")
include("HMM.jl")
include("BGMM.jl")

end # module


text:
# # Gaussian Process Classification
#
# ## Preliminary steps
#
# ### Loading necessary packages

using Plots
using HTTP, CSV
using DataFrames: DataFrame
using AugmentedGaussianProcesses

# ### Loading the banana dataset from OpenML
data = HTTP.get("https://www.openml.org/data/get_csv/1586217/phpwRjVjk")
data = CSV.read(data.body, DataFrame)
data.Class[data.Class .== 2] .= -1
data = Matrix(data)
X = data[:, 1:2]
Y = data[:, end];

# ### We create a function to visualize the data

function plot_data(X, Y; size=(300, 500))
    return Plots.scatter(
        eachcol(X)...; group=Y, alpha=0.2, markerstrokewidth=0.0, lab="", size=size
    )
end
plot_data(X, Y; size=(500, 500))

# ### Run sparse classification with increasing number of inducing points
Ms = [4, 8, 16, 32, 64]
models = Vector{AbstractGP}(undef, length(Ms) + 1)
kernel = SqExponentialKernel() ∘ ScaleTransform(1.0)
for (i, num_inducing) in enumerate(Ms)
    @info "Training with $(num_inducing) points"
    m = SVGP(
        X,
        Y,
        kernel,
        LogisticLikelihood(),
        AnalyticVI(),
        num_inducing;
        optimiser=false,
        Zoptimiser=false,
    )
    @time train!(m, 20)
    models[i] = m
end
# ### Running the full model
@info "Running full model"
mfull = VGP(X, Y, kernel, LogisticLikelihood(), AnalyticVI(); optimiser=false)
@time train!(mfull, 5)
models[end] = mfull

# ### We create a prediction and plot function on a grid
function compute_grid(model, n_grid=50)
    mins = [-3.25, -2.85]
    maxs = [3.65, 3.4]
    x_lin = range(mins[1], maxs[1]; length=n_grid)
    y_lin = range(mins[2], maxs[2]; length=n_grid)
    x_grid = Iterators.product(x_lin, y_lin)
    y_grid, _ = proba_y(model, vec(collect.(x_grid)))
    return y_grid, x_lin, y_lin
end

function plot_model(model, X, Y, title=nothing; size=(300, 500))
    n_grid = 50
    y_pred, x_lin, y_lin = compute_grid(model, n_grid)
    title = if isnothing(title)
        (model isa SVGP ? "M = $(AGP.dim(model[1]))" : "full")
    else
        title
    end
    p = plot_data(X, Y; size=size)
    Plots.contour!(
        p,
        x_lin,
        y_lin,
        reshape(y_pred, n_grid, n_grid)';
        cbar=false,
        levels=[0.5],
        fill=false,
        color=:black,
        linewidth=2.0,
        title=title,
    )
    if model isa SVGP
        Plots.scatter!(
            p, eachrow(hcat(AGP.Zview(model[1])...))...; msize=2.0, color="black", lab=""
        )
    end
    return p
end;

# ### Now run the prediction for every model and visualize the differences
Plots.plot(
    plot_model.(models, Ref(X), Ref(Y))...; layout=(1, length(models)), size=(1000, 200)
)

# ## Bayesian SVM vs Logistic
# ### We now create a model with the Bayesian SVM likelihood

mbsvm = VGP(X, Y, kernel, BayesianSVM(), AnalyticVI(); optimiser=false)
@time train!(mbsvm, 5)
# ### And compare it with the Logistic likelihood
Plots.plot(
    plot_model.(
        [models[end], mbsvm], Ref(X), Ref(Y), ["Logistic", "BSVM"]; size=(500, 500)
    )...;
    layout=(1, 2),
)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayOnDemandProbe <: SwaggerModel
    protocol::Any # spec type: Union{ Nothing, ApplicationGatewayProtocol } # spec name: protocol
    host::Any # spec type: Union{ Nothing, String } # spec name: host
    path::Any # spec type: Union{ Nothing, String } # spec name: path
    timeout::Any # spec type: Union{ Nothing, Int32 } # spec name: timeout
    pickHostNameFromBackendHttpSettings::Any # spec type: Union{ Nothing, Bool } # spec name: pickHostNameFromBackendHttpSettings
    match::Any # spec type: Union{ Nothing, ApplicationGatewayProbeHealthResponseMatch } # spec name: match
    backendAddressPool::Any # spec type: Union{ Nothing, SubResource } # spec name: backendAddressPool
    backendHttpSettings::Any # spec type: Union{ Nothing, SubResource } # spec name: backendHttpSettings

    function ApplicationGatewayOnDemandProbe(;protocol=nothing, host=nothing, path=nothing, timeout=nothing, pickHostNameFromBackendHttpSettings=nothing, match=nothing, backendAddressPool=nothing, backendHttpSettings=nothing)
        o = new()
        validate_property(ApplicationGatewayOnDemandProbe, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(ApplicationGatewayOnDemandProbe, Symbol("host"), host)
        setfield!(o, Symbol("host"), host)
        validate_property(ApplicationGatewayOnDemandProbe, Symbol("path"), path)
        setfield!(o, Symbol("path"), path)
        validate_property(ApplicationGatewayOnDemandProbe, Symbol("timeout"), timeout)
        setfield!(o, Symbol("timeout"), timeout)
        validate_property(ApplicationGatewayOnDemandProbe, Symbol("pickHostNameFromBackendHttpSettings"), pickHostNameFromBackendHttpSettings)
        setfield!(o, Symbol("pickHostNameFromBackendHttpSettings"), pickHostNameFromBackendHttpSettings)
        validate_property(ApplicationGatewayOnDemandProbe, Symbol("match"), match)
        setfield!(o, Symbol("match"), match)
        validate_property(ApplicationGatewayOnDemandProbe, Symbol("backendAddressPool"), backendAddressPool)
        setfield!(o, Symbol("backendAddressPool"), backendAddressPool)
        validate_property(ApplicationGatewayOnDemandProbe, Symbol("backendHttpSettings"), backendHttpSettings)
        setfield!(o, Symbol("backendHttpSettings"), backendHttpSettings)
        o
    end
end # type ApplicationGatewayOnDemandProbe

const _property_map_ApplicationGatewayOnDemandProbe = Dict{Symbol,Symbol}(Symbol("protocol")=>Symbol("protocol"), Symbol("host")=>Symbol("host"), Symbol("path")=>Symbol("path"), Symbol("timeout")=>Symbol("timeout"), Symbol("pickHostNameFromBackendHttpSettings")=>Symbol("pickHostNameFromBackendHttpSettings"), Symbol("match")=>Symbol("match"), Symbol("backendAddressPool")=>Symbol("backendAddressPool"), Symbol("backendHttpSettings")=>Symbol("backendHttpSettings"))
const _property_types_ApplicationGatewayOnDemandProbe = Dict{Symbol,String}(Symbol("protocol")=>"ApplicationGatewayProtocol", Symbol("host")=>"String", Symbol("path")=>"String", Symbol("timeout")=>"Int32", Symbol("pickHostNameFromBackendHttpSettings")=>"Bool", Symbol("match")=>"ApplicationGatewayProbeHealthResponseMatch", Symbol("backendAddressPool")=>"SubResource", Symbol("backendHttpSettings")=>"SubResource")
Base.propertynames(::Type{ ApplicationGatewayOnDemandProbe }) = collect(keys(_property_map_ApplicationGatewayOnDemandProbe))
Swagger.property_type(::Type{ ApplicationGatewayOnDemandProbe }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayOnDemandProbe[name]))}
Swagger.field_name(::Type{ ApplicationGatewayOnDemandProbe }, property_name::Symbol) =  _property_map_ApplicationGatewayOnDemandProbe[property_name]

function check_required(o::ApplicationGatewayOnDemandProbe)
    true
end

function validate_property(::Type{ ApplicationGatewayOnDemandProbe }, name::Symbol, val)
end


text:
module IndirectArrays

export IndirectArray

"""
    IndirectArray(index, values)

creates an array `A` where the values are looked up in the value table,
`values`, using the `index`.  Concretely, `A[i,j] =
values[index[i,j]]`.
"""
immutable IndirectArray{T,N,I<:Integer} <: AbstractArray{T,N}
    index::Array{I,N}
    values::Vector{T}

    @inline function IndirectArray(index, values)
        # The typical logic for testing bounds and then using
        # @inbounds will not check whether index is inbounds for
        # values. So we had better check this on construction.
        @boundscheck checkbounds(values, index)
        new(index, values)
    end
end
Base.@propagate_inbounds IndirectArray{T,N,I<:Integer}(index::Array{I,N},values::Vector{T}) = IndirectArray{T,N,I}(index,values)

Base.size(A::IndirectArray) = size(A.index)
Base.linearindexing(A::IndirectArray) = Base.LinearFast()

@inline function Base.getindex(A::IndirectArray, i::Int)
    @boundscheck checkbounds(A.index, i)
    @inbounds idx = A.index[i]
    @boundscheck checkbounds(A.values, idx)
    @inbounds ret = A.values[idx]
    ret
end

@inline function Base.getindex{T,N}(A::IndirectArray{T,N}, I::Vararg{Int,N})
    @boundscheck checkbounds(A.index, I...)
    @inbounds idx = A.index[I...]
    @boundscheck checkbounds(A.values, idx)
    @inbounds ret = A.values[idx]
    ret
end

end # module


text: const ASSET_FINGERPRINT = "444234794cbd4f3ef32d"

max_stars_repo_name: essenciary/genie-demo-tweet-stats

max_issues_repo_name: essenciary/genie-demo-tweet-stats

max_forks_repo_name: essenciary/genie-demo-tweet-stats


text:
__precompile__()

module Twitter

##############################################################################
##
## Dependencies
##
##############################################################################

using DataFrames, OAuth, JSON
import DataFrames: DataFrame

##############################################################################
##
## Exported methods and types
##
##############################################################################

export  get_direct_messages,
        get_direct_messages_sent,
        get_direct_messages_show,
        post_direct_messages_destroy,
        post_direct_messages_send,
        get_favorites_list,
        post_favorites_destroy,
        post_favorites_create,
        get_friendships_no_retweets,
        get_friends_ids,
        get_followers_ids,
        get_friendships_lookup,
        get_friendships_incoming,
        get_friendships_outgoing,
        post_friendships_create,
        post_friendships_destroy,
        post_friendships_update,
        get_friendships_show,
        get_friends_list,
        get_followers_list,
        twitterauth,
        get_geo_id_place_id,
        get_geo_reverse_geocode,
        get_geo_search,
        get_geo_similar_places,
        post_geo_place,
        get_help_configuration,
        get_help_languages,
        get_help_privacy,
        get_help_tos,
        get_application_rate_limit_status,
        get_lists,
        get_lists_statuses,
        post_lists_members_destroy,
        get_lists_memberships,
        get_lists_subscribers,
        post_lists_subscribers_create,
        get_lists_subscribers_show,
        post_lists_subscribers_destroy,
        post_lists_members_createall,
        get_lists_members_show,
        get_lists_members,
        post_lists_members_create,
        post_lists_destroy,
        post_lists_update,
        post_lists_create,
        get_lists_show,
        get_lists_subscriptions,
        post_lists_members_destroy_all,
        get_lists_ownerships,
        get_saved_searches_list,
        get_search_tweets,
        post_users_report_spam,
        get_user_suggestions_slug,
        get_user_suggestions,
        get_user_suggestions_slug_members,
        get_mentions_timeline,
        get_user_timeline,
        get_home_timeline,
        get_retweets_of_me,
        get_trends_place,
        get_trends_available,
        get_trends_closest,
        get_retweets_id,
        get_single_tweet_id,
        post_destroy_single_tweet_id,
        post_status_update,
        post_status_retweet_id,
        get_oembed,
        get_retweeters_id,
        get_account_settings,
        get_verify_credentials,
        post_account_settings,
        post_update_delivery_device,
        post_update_profile,
        post_update_profile_background,
        post_update_profile_image,
        get_blocks_list,
        get_blocks_ids,
        post_blocks_create,
        post_blocks_destroy,
        get_users_lookup,
        get_users_show,
        get_users_search,
        get_users_contributees,
        get_users_contributors,
        post_account_remove_profile_banner,
        post_account_update_profile_banner,
        get_profile_banner,
        DataFrame,
        Lists,
        Users,
        Places,
        Tweets

##############################################################################
##
## Load files
##
##############################################################################

include("general.jl")
include("macrogenerated.jl")
include("othermethods.jl")
include("types.jl")
include("dataframe.jl")

end # module


text:
### AMSMMALA state subtypes

## MuvAMSMMALAState holds the internal state ("local variables") of the AMSMMALA sampler for multivariate parameters

type MuvAMSMMALAState <: MuvPSMMALAState
  pstate::ParameterState{Continuous, Multivariate}
  tune::MCTunerState
  sqrttunestep::Real
  ratio::Real
  μ::RealVector
  lastmean::RealVector
  secondlastmean::RealVector
  newinvtensor::RealMatrix
  oldinvtensor::RealMatrix
  cholinvtensor::RealLowerTriangular
  newfirstterm::RealVector
  oldfirstterm::RealVector
  presentupdatetensor::Bool
  pastupdatetensor::Bool
  count::Integer
  updatetensorcount::Integer

  function MuvAMSMMALAState(
    pstate::ParameterState{Continuous, Multivariate},
    tune::MCTunerState,
    sqrttunestep::Real,
    ratio::Real,
    μ::RealVector,
    lastmean::RealVector,
    secondlastmean::RealVector,
    newinvtensor::RealMatrix,
    oldinvtensor::RealMatrix,
    cholinvtensor::RealLowerTriangular,
    newfirstterm::RealVector,
    oldfirstterm::RealVector,
    presentupdatetensor::Bool,
    pastupdatetensor::Bool,
    count::Integer,
    updatetensorcount::Integer
  )
    if !isnan(ratio)
      @assert 0 < ratio < 1 "Acceptance ratio should be between 0 and 1"
      @assert sqrttunestep > 0 "Square root of tuned drift step is not positive"
    end
    new(
      pstate,
      tune,
      sqrttunestep,
      ratio,
      μ,
      lastmean,
      secondlastmean,
      newinvtensor,
      oldinvtensor,
      cholinvtensor,
      newfirstterm,
      oldfirstterm,
      presentupdatetensor,
      pastupdatetensor,
      count,
      updatetensorcount
    )
  end
end

MuvAMSMMALAState(pstate::ParameterState{Continuous, Multivariate}, tune::MCTunerState=PSMMALAMCTune()) =
  MuvAMSMMALAState(
  pstate,
  tune,
  NaN,
  NaN,
  Array(eltype(pstate), pstate.size),
  Array(eltype(pstate), pstate.size),
  Array(eltype(pstate), pstate.size),
  Array(eltype(pstate), pstate.size, pstate.size),
  Array(eltype(pstate), pstate.size, pstate.size),
  RealLowerTriangular(Array(eltype(pstate), pstate.size, pstate.size)),
  Array(eltype(pstate), pstate.size),
  Array(eltype(pstate), pstate.size),
  false,
  false,
  0,
  0
)

### Metropolis-adjusted Langevin Algorithm (AMSMMALA)

immutable AMSMMALA <: PSMMALA
  driftstep::Real
  update!::Function
  transform::Union{Function, Void}
  t0::Integer

  function AMSMMALA(driftstep::Real, update!::Function, transform::Union{Function, Void}, t0::Integer)
    @assert driftstep > 0 "Drift step is not positive"
    @assert t0 > 0 "t0 is not positive"
    new(driftstep, update!, transform, t0)
  end
end

AMSMMALA(driftstep::Real=1.; update::Function=rand_update!, transform::Union{Function, Void}=nothing, t0::Integer=3) =
  AMSMMALA(driftstep, update, transform, t0)

### Initialize AMSMMALA sampler

## Initialize parameter state

function initialize!(
  pstate::ParameterState{Continuous, Multivariate},
  parameter::Parameter{Continuous, Multivariate},
  sampler::AMSMMALA,
  outopts::Dict
)
  parameter.uptotensorlogtarget!(pstate)
  if sampler.transform != nothing
    pstate.tensorlogtarget = sampler.transform(pstate.tensorlogtarget)
  end
  pstate.diagnosticvalues[1] = false
  @assert isfinite(pstate.logtarget) "Log-target not finite: initial value out of support"
  @assert all(isfinite(pstate.gradlogtarget)) "Gradient of log-target not finite: initial values out of support"
  @assert all(isfinite(pstate.tensorlogtarget)) "Tensor of log-target not finite: initial values out of support"

  if !isempty(outopts[:diagnostics])
    pstate.diagnostickeys = copy(outopts[:diagnostics])
    pstate.diagnosticvalues = Array(Any, length(pstate.diagnostickeys))
  end
end

## Initialize AMSMMALA state

tuner_state(sampler::AMSMMALA, tuner::PSMMALAMCTuner) =
  PSMMALAMCTune(
    BasicMCTune(NaN, 0, 0, tuner.smmalatuner.period),
    BasicMCTune(NaN, 0, 0, tuner.malatuner.period),
    BasicMCTune(sampler.driftstep, 0, 0, tuner.totaltuner.period)
  )

function sampler_state(
  parameter::Parameter{Continuous, Multivariate},
  sampler::AMSMMALA,
  tuner::MCTuner,
  pstate::ParameterState{Continuous, Multivariate},
  vstate::VariableStateVector
)
  sstate = MuvAMSMMALAState(generate_empty(pstate, parameter.diffmethods, parameter.diffopts), tuner_state(sampler, tuner))
  sstate.sqrttunestep = sqrt(sampler.driftstep)
  sstate.lastmean = copy(pstate.value)
  sstate.oldinvtensor = inv(pstate.tensorlogtarget)
  sstate.cholinvtensor = ctranspose(chol(Hermitian(sstate.oldinvtensor)))
  sstate.oldfirstterm = sstate.oldinvtensor*pstate.gradlogtarget
  sstate.presentupdatetensor = true
  sstate.pastupdatetensor = false
  sstate
end

### Reset AMSMMALA sampler

## Reset parameter state

function reset!(
  pstate::ParameterState{Continuous, Multivariate},
  x::RealVector,
  parameter::Parameter{Continuous, Multivariate},
  sampler::AMSMMALA
)
  pstate.value = copy(x)
  parameter.uptotensorlogtarget!(pstate)
  pstate.diagnosticvalues[1] = false
end

## Reset sampler state

function reset!(
  sstate::MuvAMSMMALAState,
  pstate::ParameterState{Continuous, Multivariate},
  parameter::Parameter{Continuous, Multivariate},
  sampler::AMSMMALA,
  tuner::MCTuner
)
  reset!(sstate.tune, sampler, tuner)
  sstate.sqrttunestep = sqrt(sampler.driftstep)
  sstate.lastmean = copy(pstate.value)
  sstate.oldinvtensor = inv(pstate.tensorlogtarget)
  sstate.cholinvtensor = chol(sstate.oldinvtensor, Val{:L})
  sstate.oldfirstterm = sstate.oldinvtensor*pstate.gradlogtarget
  sstate.presentupdatetensor = true
  sstate.pastupdatetensor = false
end

Base.show(io::IO, sampler::AMSMMALA) = print(io, "AMSMMALA sampler: drift step = $(sampler.driftstep)")

Base.show(io::IO, ::MIME"text/plain", sampler::AMSMMALA) = show(io, sampler)


text:
using Mustache
using GraphLayout
using DataFrames
using Colors

function Base.writemime(f::IO, ::MIME"image/svg+xml", nw::CommunityNetwork)

    A = adjacency_matrix(nw.graph)
    nx, ny = layout_spring_adj(full(A))

    scale = sqrt(length(nw.graph.vertices))/sqrt(200)
    trans(x) = scale*250.0*x+250.0

    #=
    first = Int[l.first for l in nw.graph.edges]
    second = Int[l.second for l in nw.graph.edges]

    edge_df=DataFrame(x1s=trans(nx[first]),
                    x2s=trans(nx[second]),
                    y1s=trans(ny[first]),
                    y2s=trans(ny[second]))

    node_df=DataFrame(xs=trans(nx),
                      ys=trans(ny))
    =#

    gedges = collect(nw.graph.edges)
    edges = Array(Dict{Any,Any}, length(gedges))
    for i in 1:length(gedges)
        edges[i] = {:x1s=>trans(nx[gedges[i].first]),
                    :x2s=>trans(nx[gedges[i].second]),
                    :y1s=>trans(ny[gedges[i].first]),
                    :y2s=>trans(ny[gedges[i].second])}
    end

    nodes = Array(Dict{Any,Any}, length(nw.graph.vertices))
    for i in 1:length(nw.graph.vertices)
        nodes[i] = {:xs=>trans(nx[i]),
                    :ys=>trans(ny[i]),
                    :color=>"hsl($(360/nw.n_comms*nw.labels[i]),100%,50%)"}
    end

    template = readall(joinpath(dirname(@__FILE__()), "graph_template.svg"))
    #XXX tx and ty are a hack because I always used node 1 as the target
    svg_string = render(template, edges=edges, nodes=nodes, tx=trans(nx[1]), ty=trans(ny[1]))

    # # if you don't want to reload the svg every time:
    # filename = joinpath(dirname(@__FILE__()), "graph_template.svg")
    # svg_string = render_from_file(filename, edges=edges, nodes=nodes)

    # debugging
    # println(svg_string)

    println(f,svg_string)
end


text:
using Test
using FTCTests

@testset "run_sim" begin
    dir_log = "__data__"
    dir_figure = "__figures__"
    case_number = 1
    mkpath(dir_log)
    method = :adaptive
    args_multicopter = ()  # default
    multicopter = LeeHexacopter()
    fault = FaultSet(
                     LoE(3.0, 1, 0.9),
                    )
    fdi = DelayFDI(0.1)  # 0.1 sec delay
    t0, tf = 0.0, 0.1
    θs = [[0, 0, 0.0]]  # constant position tracking
    traj_des = Bezier(θs, t0, tf)
    run_sim(method, args_multicopter, multicopter, fault, fdi, traj_des, dir_log, case_number;
            t0=t0, tf=tf, savestep=0.01, dir_figure=dir_figure)
end


text:
#!/usr/bin/env julia

const DATA_DIR   = "../18_data"
const DEBUG      = false
const INPUT_FILE = "input.txt"

const OP_PLUS  = '+'
const OP_TIMES = '*'
const OPS = [OP_PLUS, OP_TIMES]

const PAREN_LEFT  = '('
const PAREN_RIGHT = ')'
const PARENS      = [PAREN_LEFT, PAREN_RIGHT]

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

# Note: Form parenthetical groups from op eval
function control_precedence(c)
    if c == PAREN_LEFT
        return 0
    elseif c == OP_PLUS
        return 2
    elseif c == OP_TIMES
        return 1
    else
        error("Unrecognized operator")
    end
end

function op_eval(a, op, b)
    if op == '+'
        return a + b
    elseif op == '*'
        return a * b
    end
end

function dprintln(s)
   if DEBUG
       println(s)
   end
end

function top_eval!(vals, controls)
    val2 = pop!(vals)
    val1 = pop!(vals)
    op = pop!(controls)
    opval = op_eval(val1, op, val2)
    push!(vals, opval)
    dprintln("\tINFO: top_eval!: val1=$(val1), op=$(op), val2=$(val2), opval=$(opval)")
end

function expr_eval(expr::String)
    dprintln("INFO: expr_eval: Evaluating $(expr)")
    vals = []
    controls = []  # Ops & lparens, but not rparens
    loc = 1
    while loc <= length(expr)
        c = expr[loc]
        if isspace(c)
            loc += 1
            continue
        end
        dprintln("INFO:    Processing token $(c)")

        if isdigit(c)
            val = parse(Int, c)
            ### loc += 1
            ### while loc <= length(expr) && isdigit(expr[loc])
            ###     val = (val * 10) + parse(Int, expr[loc])
            ###     loc += 1
            ### end
            ### loc -= 1
            push!(vals, val)
            dprintln("\tINFO: LOOP(digit): Pushed val=$(val)")
        elseif c == PAREN_LEFT
            push!(controls, c)
            dprintln("\tINFO: LOOP(lparen): Pushed control=$(c)")
        elseif c == PAREN_RIGHT  # Eval everything back to last lparen
            while !isempty(controls) && last(controls) != PAREN_LEFT
                top_eval!(vals, controls)
            end
            pop!(controls) # Remove last lparen
            dprintln("\tINFO: LOOP(rparen): Evaluated backlog and popped rparen")
        else
            @assert c in OPS

            # Evaluate previous higher-priority ops before this one
            while !isempty(controls) && control_precedence(last(controls)) >= control_precedence(c)
                # println("\tIn ops while loop: controls=$(controls)")
                top_eval!(vals, controls)
            end
            push!(controls, c)
            dprintln("\tINFO: LOOP(op): Evaluated lower-precedence ops and pushed $(c)")
        end
        loc += 1
    end

    dprintln("\tINFO: Performing clean-up")
    while !isempty(controls)
        top_eval!(vals, controls)
    end

    @assert length(vals) == 1
    @assert isempty(controls)
    return last(vals)
end

function test_eval(expr::String)
    println("Expr: $(expr)")
    val = expr_eval(expr)
    println("\tComputed value: $(val)")
    println("="^20)
end

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

function main(filename::String)
    lines::Array{String} = readlines(filename)
    total = 0
    for line in lines
        dprintln("="^20)
        val = expr_eval(line)
        dprintln("val=$(val)")
        total += val
    end
    println("Total of all exprs=$(total)")
end

main("$DATA_DIR/$INPUT_FILE")

max_stars_repo_name: jaycoskey/AdventOfCode

max_issues_repo_name: jaycoskey/AdventOfCode

max_forks_repo_name: jaycoskey/AdventOfCode


text:
## UC Model Ref
branches = Dict{Symbol, DeviceModel}()
services = Dict{Symbol, ServiceModel}()
devices = Dict(
    :Generators => DeviceModel(ThermalStandard, ThermalBasicUnitCommitment),
    :Ren => DeviceModel(RenewableDispatch, FixedOutput),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    :ILoads => DeviceModel(InterruptibleLoad, StaticPowerLoad),
)
template_basic_uc =
    OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)

branches = Dict()
services = Dict()
devices = Dict(
    :Generators => DeviceModel(ThermalStandard, ThermalStandardUnitCommitment),
    :Ren => DeviceModel(RenewableDispatch, FixedOutput),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    :ILoads => DeviceModel(InterruptibleLoad, StaticPowerLoad),
)
template_standard_uc =
    OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)

## ED Model Ref
branches = Dict()
services = Dict()
devices = Dict(
    :Generators => DeviceModel(ThermalStandard, ThermalDispatchNoMin),
    :Ren => DeviceModel(RenewableDispatch, RenewableFullDispatch),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    :ILoads => DeviceModel(InterruptibleLoad, DispatchablePowerLoad),
)
template_ed = OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)
template_ed_ptdf = OperationsProblemTemplate(StandardPTDFModel, devices, branches, services)

## UC with services Model Ref
branches = Dict()
services = Dict(
    :ReserveUp => ServiceModel(VariableReserve{ReserveUp}, RangeReserve),
    :ReserveDown => ServiceModel(VariableReserve{ReserveDown}, RangeReserve),
)
devices = Dict(
    :Generators => DeviceModel(ThermalStandard, ThermalBasicUnitCommitment),
    :Ren => DeviceModel(RenewableDispatch, FixedOutput),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    :ILoads => DeviceModel(InterruptibleLoad, StaticPowerLoad),
)
template_basic_uc_svc =
    OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)

# UC with Hydro Model Ref
branches = Dict()
services = Dict()
devices = Dict(
    :Generators => DeviceModel(ThermalStandard, ThermalBasicUnitCommitment),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    :HydroEnergyReservoir => DeviceModel(HydroEnergyReservoir, HydroDispatchRunOfRiver),
)
template_hydro_basic_uc =
    OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)

branches = Dict()
services = Dict()
devices = Dict(
    :Generators => DeviceModel(ThermalStandard, ThermalStandardUnitCommitment),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    :HydroEnergyReservoir => DeviceModel(HydroEnergyReservoir, HydroDispatchRunOfRiver),
)
template_hydro_standard_uc =
    OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)

## ED with Hydro Model Ref
branches = Dict()
services = Dict()
devices = Dict(
    :Generators => DeviceModel(ThermalStandard, ThermalDispatchNoMin),
    :Ren => DeviceModel(RenewableDispatch, RenewableFullDispatch),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    :ILoads => DeviceModel(InterruptibleLoad, DispatchablePowerLoad),
    :HydroEnergyReservoir =>
        DeviceModel(HydroEnergyReservoir, HydroDispatchReservoirBudget),
)
template_hydro_ed =
    OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)

# UC with Hydro Model Ref
branches = Dict()
services = Dict()
devices = Dict(
    :Generators => DeviceModel(ThermalStandard, ThermalStandardUnitCommitment),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    :HydroEnergyReservoir => DeviceModel(HydroEnergyReservoir, HydroDispatchRunOfRiver),
)
template_pwl_standard_uc =
    OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)

## ED with Hydro Model Ref
branches = Dict()
services = Dict()
devices = Dict(
    :Generators => DeviceModel(ThermalStandard, ThermalRampLimited),
    :Ren => DeviceModel(RenewableDispatch, RenewableFullDispatch),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    :ILoads => DeviceModel(InterruptibleLoad, DispatchablePowerLoad),
    :HydroEnergyReservoir =>
        DeviceModel(HydroEnergyReservoir, HydroDispatchReservoirBudget),
)
template_pwl_ed =
    OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)

function PSI._jump_value(int::Int)
    @warn("This is for testing purposes only.")
    return int
end

function _test_plain_print_methods(list::Array)
    for object in list
        normal = repr(object)
        io = IOBuffer()
        show(io, "text/plain", object)
        grabbed = String(take!(io))
        @test !isnothing(grabbed)
    end
end

function _test_html_print_methods(list::Array)
    for object in list
        normal = repr(object)
        io = IOBuffer()
        show(io, "text/html", object)
        grabbed = String(take!(io))
        @test !isnothing(grabbed)
    end
end

struct FakeStagesStruct
    stages::Dict{Int, Int}
end
function Base.show(io::IO, struct_stages::FakeStagesStruct)
    PSI._print_inter_stages(io, struct_stages.stages)
    println(io, "\n\n")
    PSI._print_intra_stages(io, struct_stages.stages)
end

branches = Dict()
services = Dict()
devices = Dict(
    :Generators => DeviceModel(ThermalStandard, ThermalStandardUnitCommitment),
    :Ren => DeviceModel(RenewableDispatch, RenewableFullDispatch),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    :ILoads => DeviceModel(InterruptibleLoad, DispatchablePowerLoad),
    :HydroEnergyReservoir =>
        DeviceModel(HydroEnergyReservoir, HydroDispatchReservoirStorage),
)
template_hydro_st_uc =
    OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)

## ED with HydroEnergyReservoir Model Ref
branches = Dict()
services = Dict()
devices = Dict(
    :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
    :Ren => DeviceModel(RenewableDispatch, RenewableFullDispatch),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    :ILoads => DeviceModel(InterruptibleLoad, DispatchablePowerLoad),
    :HydroEnergyReservoir =>
        DeviceModel(HydroEnergyReservoir, HydroDispatchReservoirStorage),
)
template_hydro_st_ed =
    OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)
#=
## UC Model Ref
branches = Dict(:L => DeviceModel(Line, StaticLine),
                                     :T => DeviceModel(Transformer2W, StaticTransformer),
                                     :TT => DeviceModel(TapTransformer, StaticTransformer),
                                     :dc_line => DeviceModel(HVDCLine, HVDCDispatch))

services = Dict()

devices = Dict(:Generators => DeviceModel(ThermalStandard, ThermalStandardUnitCommitment),
                                    :Ren => DeviceModel(RenewableDispatch, RenewableFullDispatch),
                                    :Loads =>  DeviceModel(PowerLoad, StaticPowerLoad),
                                    :ILoads =>  DeviceModel(InterruptibleLoad, StaticPowerLoad))

template_basic_uc= OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)

## ED Model Ref
branches = Dict(:L => DeviceModel(Line, StaticLine),
                                     :T => DeviceModel(Transformer2W, StaticTransformer),
                                     :TT => DeviceModel(TapTransformer, StaticTransformer),
                                     :dc_line => DeviceModel(HVDCLine, HVDCDispatch))

services = Dict()

devices = Dict(:Generators => DeviceModel(ThermalStandard, ThermalDispatch, SemiContinuousFF(PSI.ACTIVE_POWER, PSI.ON),
                                    :Ren => DeviceModel(RenewableDispatch, RenewableFullDispatch),
                                    :Loads =>  DeviceModel(PowerLoad, StaticPowerLoad),
                                    :ILoads =>  DeviceModel(InterruptibleLoad, InterruptiblePowerLoad,))

template_ed= OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)
=#


text:
module FindingPi

export Point, is_inside, find_pi, random_point, generate_points

type Point
    x::Real
    y::Real
end

# r^2 = x^2 + y^2
is_inside(x, y, radius) = x^2 + y^2 <= radius^2
is_inside(point, radius) = is_inside(point.x, point.y, radius)
is_inside(point) = is_inside(point, 1)

# area_circle = pi*r^2
# area_square = (2r)^2
#  => c / s = pi / 4
find_pi(inside, total) = 4 * inside / total

function random_point()
    x = (rand() * 2) - 1
    y = (rand() * 2) - 1
    Point(x, y)
end

generate_points(n) = [random_point() for i=1:n]

end

module TestingPi

using FindingPi
using Base.Test

# is_inside
@test is_inside(0, 0, 1)
@test !is_inside(1, 1, 1)
@test is_inside(1/sqrt(2), 1/sqrt(2), 1)
@test is_inside(Point(0,0))

# find_pi
@test find_pi(100, 100) == 4

# generate_points
@test typeof(random_point()) == Point
@test length(generate_points(4)) == 4

end

module FindingPiSim

using FindingPi

function trial(n)
    points = generate_points(n)
    points_inside = map(is_inside, points)
    inside_count = sum(points_inside)
    pi_estimate = find_pi(inside_count, n)
end

function batch(range, iterations)
    results = Array(Real, (length(range) * iterations, 2))
    i = 1
    for n = range
        for j = 1:iterations
            results[i, :] = [n, trial(n)]
            i += 1
        end
    end
    results
end

end

# Analysis with Gadfly and DataFrames
# -----------------------------------

FindingPiSim.trial(1000000)

Pkg.add("Gadfly")
Pkg.add("DataFrames")

using Gadfly
using DataFrames

range = 100:10:1000
iterations = 5

total_points = sum(range) * iterations
println("$total_points points")
t1 = time()

res = FindingPiSim.batch(range, iterations)

t2 = time()
println("total time: $(t2 - t1)")
println("time per point: $((t2 - t1) / total_points)")

df = DataFrame(n=res[:, 1], pi_est=res[:, 2])
df[:error] = pi - df[:pi_est]

# A simple plot:
plot(df, x="n", y="pi_est", color="error", Geom.point)

# A more complex plot, with a theme for point size,
# and a layer to mark the true value of pi with a line:
#
# t = Theme(
#     default_point_size=Compose.Measure(0.7),
#     highlight_width=Compose.Measure(0.1)
# )
#
# plot(
#     layer(x->3.1415, range[1], range[end]),
#     layer(df, x="n", y="pi_est", color="error", Geom.point, t)
# )

max_stars_repo_name: swans-one/IntroToJulia

max_issues_repo_name: swans-one/IntroToJulia

max_forks_repo_name: swans-one/IntroToJulia


text:
export create_grid

"""
    create_grid(bounds[, steps])

Create a Cartesian grid by discretizing the line segment, rectangle, or box
defined by `bounds` in increments of `steps`.  If `steps` is not defined, a
default step size of 1.0 will be used in all directions.

# Arguments
- `bounds::Array{<:Number}`: Array that defines the min/max values of a line
  segment (1D), rectangle (2D), or box (3D) that the grid will be placed in.
  The first row of the array defines the minimum values, the second row defines
  the maximum values.
- `steps::Array{Union{AbstractFloat, Integer}}`: Array that defines how many
  steps in each direction.  If the array values in `steps` are of type
  `<:Integer` integers then those values define how many points are used in
  each direction.  If the array values in `steps` are of type `<:AbstractFloat`
  numbers then those values define the step size in each direction.

# Examples
```jldoctest
julia> create_grid([0; 2], [2])
2-element Array{Float64,1}:
 0.0
 2.0
```

```jldoctest
julia> create_grid([0 0; 2 2], [2.0 1.0])
6×2 Array{Float64,2}:
 0.0  0.0
 0.0  1.0
 0.0  2.0
 2.0  0.0
 2.0  1.0
 2.0  2.0
```

```jldoctest
julia> create_grid([0 0 0; 2 2 2], [2 2 2])
8×3 Array{Float64,2}:
 0.0  0.0  0.0
 0.0  0.0  2.0
 0.0  2.0  0.0
 0.0  2.0  2.0
 2.0  0.0  0.0
 2.0  0.0  2.0
 2.0  2.0  0.0
 2.0  2.0  2.0
```
"""
function create_grid(bounds::Array{<:Number}, steps::Array{T}) where {T<:Union{AbstractFloat, Integer}}

    if length(size(bounds)) == 1
        dimension = 1
    else
        dimension = size(bounds)[2]
    end

    if dimension ≠ length(steps)
        throw(ArgumentError("""Dimension mismatch in arguments.  bounds indicates
                dimension = $dimension. steps indicates dimension = $(length(steps))"""))
    end

    if T <: AbstractFloat
        my_range = (start, step, stop) -> start:step:stop
    elseif T <: Integer
        my_range = (start, step, stop) -> LinRange(start, stop, step)
    end

    if dimension == 1
        xrange = my_range(bounds[1], steps[1], bounds[2])
        number_of_points = length(xrange)
        x = [ i for i in xrange ]
    elseif dimension == 2
        xrange = my_range(bounds[1,1], steps[1], bounds[2, 1])
        yrange = my_range(bounds[1,2], steps[2], bounds[2, 2])
        x = [ x for x in xrange for j in 1:length(yrange) ]
        y = [ y for i in 1:length(xrange) for y in yrange ]
        [x y]
    elseif dimension == 3
        xrange = my_range(bounds[1,1], steps[1], bounds[2, 1])
        yrange = my_range(bounds[1,2], steps[2], bounds[2, 2])
        zrange = my_range(bounds[1,2], steps[3], bounds[2, 3])
        x = [ x for x in xrange for j in 1:length(yrange) for k in 1:length(zrange) ]
        y = [ y for i in 1:length(xrange) for y in yrange for k in 1:length(zrange) ]
        z = [ z for i in 1:length(xrange) for j in 1:length(yrange) for z in zrange ]
        [x y z]
    else
        throw(ArgumentError("Can't have more than 3 spatial dimensions for grid."))
    end

end

function create_grid(bounds::Array{<:Number})

    if length(size(bounds)) == 1
        dimension = 1
    else
        dimension = size(bounds)[2]
    end

    if dimension == 1
        steps = [1.0]
    elseif dimension == 2
        steps = [1.0 1.0]
    elseif dimension == 3
        steps = [1.0 1.0 1.0]
    else
        throw(ArgumentError("Can't have more than 3 spatial dimensions for grid."))
    end

    create_grid(bounds, steps)

end


text:
# Model
include(joinpath(pwd(), "models/quadruped.jl"))
model = Quadruped{Discrete, FixedTime}(n, m, d,
			  g, 0.1,
			  l_torso, d_torso, m_torso, J_torso,
			  l_thigh, d_thigh, m_thigh, J_thigh,
			  l_leg, d_leg, m_leg, J_leg,
			  l_thigh, d_thigh, m_thigh, J_thigh,
			  l_leg, d_leg, m_leg, J_leg,
			  l_thigh, d_thigh, m_thigh, J_thigh,
			  l_leg, d_leg, m_leg, J_leg,
			  l_thigh, d_thigh, m_thigh, J_thigh,
			  l_leg, d_leg, m_leg, J_leg,
			  qL, qU,
			  uL, uU,
			  nq,
			  nu,
			  nc,
			  nf,
			  nb,
			  ns,
			  idx_u,
			  idx_λ,
			  idx_b,
			  idx_ψ,
			  idx_η,
			  idx_s,
			  joint_friction)
model_ft = free_time_model(model)

# Horizon
T = 41

# Time step
tf = 0.5
h = tf / (T - 1)

# Visualize
# - Pkg.add any external deps from visualize.jl
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)

# Configurations
θ = pi / 2.5
q1 = initial_configuration(model_ft,  θ)
q1_high = initial_configuration(model_ft,  π / 3.5)

qM = copy(q1)
qM[1] -= 0.5 * model_ft.l_torso
qM[2] += 0.75
# qM[3] += pi
# qM[4] += pi
# qM[5] += pi
# qM[6] += pi
# qM[7] += pi
# qM[8] += pi
# qM[9] += pi
# qM[10] += pi
# qM[11] += pi
qT = copy(q1_high)
qT[1] -= 2 * model_ft.l_torso
# qT[3] += 2.0 * pi
# qT[4] += 2.0 * pi
# qT[5] += 2.0 * pi
# qT[6] += 2.0 * pi
# qT[7] += 2.0 * pi
# qT[8] += 2.0 * pi
# qT[9] += 2.0 * pi
# qT[10] += 2.0 * pi
# qT[11] += 2.0 * pi

visualize!(vis, model_ft, [q1])
visualize!(vis, model_ft, [qM])
visualize!(vis, model_ft, [qT])

q_ref = [linear_interpolation(q1, qM, 21)...,
    linear_interpolation(qM, qT, 21)...]

visualize!(vis, model_ft, q_ref)

# Bounds

# control
# u = (τ1..4, λ1..2, β1..4, ψ1..2, η1...4, s1)
# ul <= u <= uu
_uu = Inf * ones(model_ft.m)
_uu[model_ft.idx_u] .= Inf
_uu[end] = 5.0 * h

_ul = zeros(model_ft.m)
_ul[model_ft.idx_u] .= -Inf
_ul[end] = 0.5 * h
ul, uu = control_bounds(model_ft, T, _ul, _uu)

xl, xu = state_bounds(model_ft, T,
    x1 = [q1; q1],
    xT = [Inf * ones(nq); qT])

# Objective
include_objective(["velocity", "nonlinear_stage"])

x0 = configuration_to_state(q_ref)

# penalty on slack variable
obj_penalty = PenaltyObjective(1.0e4, model_ft.m - 1)

# quadratic tracking objective
# Σ (x - xref)' Q (x - x_ref) + (u - u_ref)' R (u - u_ref)
obj_control = quadratic_time_tracking_objective(
    [Diagonal(0.0 * ones(2 * model.nq)) for t = 1:T],
    [Diagonal([0.0 * ones(model_ft.nu)..., 0.0 * ones(model_ft.nc + model_ft.nb)..., zeros(model_ft.m - model_ft.nu - model_ft.nc - model_ft.nb)...]) for t = 1:T-1],
    [deepcopy(x0[t]) for t = 1:T],
    [zeros(model_ft.m) for t = 1:T],
    1.0)

# quadratic velocity penalty
# Σ v' Q v
obj_velocity = velocity_objective(
    [Diagonal(1.0 * ones(model_ft.nq)) for t = 1:T-1],
    model_ft.nq,
    h = h,
    idx_angle = collect([3, 4, 5, 6, 7, 8, 9, 10, 11]))

t_h = kinematics_1(model_ft, qM, body = :torso, mode = :com)[2]
θ_range = range(0.5 * π, stop = 0.5 * π + 2 * π, length = T)
function l_stage(x, u, t)
	J = 0.0

	q = view(x, 1:11)

	J += 10000.0 * (kinematics_1(model_ft, q, body = :torso, mode = :com)[2] - kinematics_1(model_ft, q_ref[t], body = :torso, mode = :com)[2])^2.0
	# J += 10000.0 * (q[3] - θ_range[t])^2.0

	return J
end

l_stage(x) = l_stage(x, nothing, T)
obj_stage = nonlinear_stage_objective(l_stage, l_stage)

obj = MultiObjective([obj_penalty,
                      obj_control,
                      obj_velocity,
                      obj_stage])

# Constraints
include_constraints(["stage", "contact", "free_time"])
con_contact = contact_constraints(model_ft, T)
con_free_time = free_time_constraints(T)

pr1_1 = kinematics_2(model, q1, body = :calf_1, mode = :ee)
pr2_1 = kinematics_2(model, q1, body = :calf_2, mode = :ee)

pf1_1 = kinematics_3(model, q1, body = :calf_3, mode = :ee)
pf2_1 = kinematics_3(model, q1, body = :calf_4, mode = :ee)

pr1_T = kinematics_2(model, qT, body = :calf_1, mode = :ee)
pr2_T = kinematics_2(model, qT, body = :calf_2, mode = :ee)

pf1_T = kinematics_3(model, qT, body = :calf_3, mode = :ee)
pf2_T = kinematics_3(model, qT, body = :calf_4, mode = :ee)

function pinned_1!(c, x, u, t)
    q = view(x, 1:11)
    c[1:2] = pr1_1 - kinematics_2(model, q, body = :calf_1, mode = :ee)
    c[3:4] = pf1_1 - kinematics_3(model, q, body = :calf_3, mode = :ee)
	c[5:6] = pr2_1 - kinematics_2(model, q, body = :calf_2, mode = :ee)
    c[7:8] = pf2_1 - kinematics_3(model, q, body = :calf_4, mode = :ee)
    nothing
end

function pinned_T!(c, x, u, t)
    q = view(x, 11 .+ (1:11))
    c[1:2] = pr1_T - kinematics_2(model, q, body = :calf_1, mode = :ee)
    c[3:4] = pf1_T - kinematics_3(model, q, body = :calf_3, mode = :ee)
	c[5:6] = pr2_T - kinematics_2(model, q, body = :calf_2, mode = :ee)
    c[7:8] = pf2_T - kinematics_3(model, q, body = :calf_4, mode = :ee)
    nothing
end

n_stage = 8
t_idx1 = vcat([t for t = 1:3]...)
t_idxT = vcat([t for t = T-2:T]...)
con_pinned1 = stage_constraints(pinned_1!, n_stage, (1:0), t_idx1)
con_pinnedT = stage_constraints(pinned_T!, n_stage, (1:0), t_idxT)

con = multiple_constraints([con_contact, con_free_time, con_pinned1, con_pinnedT])

# Problem
prob = trajectory_optimization_problem(model_ft,
               obj,
               T,
               # h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# trajectory initialization
u0 = [[1.0e-3 * rand(model_ft.m - 1); h] for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)
# z0 .+= 1.0e-3 * randn(prob.num_var)

# Solve
# include_snopt()

@time z̄ , info = solve(prob, copy(z0),
    # nlp = :SNOPT7,
    tol = 1.0e-2, c_tol = 1.0e-2,
	max_iter = 2000,
    time_limit = 60 * 3, mapl = 5)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
tf, t, h̄ = get_time(ū)

# @save joinpath(pwd(), "examples/trajectories/quadruped_backflip.jld2") x̄ ū h̄

# Visualize
visualize!(vis, model_ft, state_to_configuration(x̄), Δt = h̄[1])

visualize!(vis, model_ft,
	[[state_to_configuration(x̄)[1] for i = 1:10]..., state_to_configuration(x̄)..., [state_to_configuration(x̄)[end] for i = 1:10]...],
	Δt = ū[1][end])

# using Plots
# plot(t[1:end-1], hcat(ū...)[model_ft.idx_u, :]', linetype = :steppost,
# 	width = 2.0, label = "", xlabel= "time (s)", ylabel = "control")

# save trajectory
qm = state_to_configuration(x̄)
um = [u[model.idx_u] for u in ū]
γm = [u[model.idx_λ] for u in ū]
bm = [u[model.idx_b] for u in ū]
ψm = [u[model.idx_ψ] for u in ū]
ηm = [u[model.idx_η] for u in ū]
hm = mean(h̄)
μm = model_ft.μ
@save joinpath(@__DIR__, "quadruped_jump_v2.jld2") qm um γm bm ψm ηm μm hm

using Plots
plot(hcat(γm...)', linetype = :steppost)
plot(hcat(bm...)', linetype = :steppost)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
"""

The exported names are
$(EXPORTS)
"""
module ODEs

using DocStringExtensions

include("ODETools/ODETools.jl")

include("TransientFETools/TransientFETools.jl")

# include("DiffEqsWrappers/DiffEqsWrappers.jl")

end #module

const GridapODEs = ODEs


text:
using Test, Unicode
Base.include(@__MODULE__, joinpath(Sys.BINDIR, "..", "share", "julia", "test", "testhelpers", "FakePTYs.jl"))
import .FakePTYs: open_fake_pty

slave, master = open_fake_pty()

CTRL_C = '\x03'

# Script that we want the REPL to execute, here simply st for Pkg REPLMode
test_script = """
using ReplMaker

function parse_to_expr(s)
    quote Meta.parse(\$s) end
end

initrepl(parse_to_expr,
         prompt_text="Expr> ",
         prompt_color = :blue,
         start_key=')',
         mode_name="Expr_mode")

) x + 1

"""*CTRL_C

function run_test()
    slave, master = open_fake_pty()
    # Start a julia process
    p = run(`$(Base.julia_cmd()) --history-file=no --startup-file=no`, slave, slave, slave; wait=false)

    # Read until the prompt
    readuntil(master, "julia>", keep=true)
    done = false
    repl_output_buffer = IOBuffer()

    # A task that just keeps reading the output
    @async begin
        while true
            done && break
            write(repl_output_buffer, readavailable(master))
        end
    end

    # Execute our "script"
    for l in split(test_script, '\n'; keepempty=false)
        write(master, l, '\n')
    end

    # Let the REPL exit
    write(master, "exit()\n")
    wait(p)
    done = true

    # Gather the output
    repl_output = String(take!(repl_output_buffer))
    println(repl_output)
    return split(repl_output, '\n'; keepempty=false)
end

out = run_test();

@test out[end-5] == "\e[?2004h\r\e[0K\e[34m\e[1mExpr> \e[0m\e[0m\r\e[6C\r\e[6C\r\e[0K\e[34m\e[1mExpr> \e[0m\e[0m\r\e[6C\r\e[6C^C\r"


text:
"""
Mass transfer coefficients.

* [`diffus_air`](@ref)
* [`diffus_water`](@ref)
* [`diffus_soil_air`](@ref)
* [`diffus_soil_water`](@ref)
* [`diffus_soil`](@ref)

# References

## Diffusivity in air

* [M98] Massman, W. J. (1998). A review of the molecular diffusivities of
  H2O, CO2, CH4, CO, O3, SO2, NH3, N2O, NO, and NO2 in air, O2 and N2 near STP.
  *Atmos. Environ.*, 32(6), 1111--1127.
  <https://doi.org/10.1016/S1352-2310(97)00391-9>
* [S10] Seibt, U. et al. (2010). A kinetic analysis of leaf uptake of COS
  and its relation to transpiration, photosynthesis and carbon isotope
  fractionation. *Biogeosci.*, 7, 333--341.
  <https://doi.org/10.5194/bg-7-333-2010>

## Diffusivity in water

* [J87] Jähne, B. et al. (1987). Measurement of the diffusion
  coefficients of sparingly soluble gases in water. *J. Geophys. Res.*,
  92(C10), 10767--10776. <https://doi.org/10.1029/JC092iC10p10767>
* [WH68] Wise, D. L. and Houghton, G. (1968). Diffusion coefficients of
  neon, krypton, xenon, carbon monoxide and nitric oxide in water at 10-6O C.
  *Chem. Eng. Sci.*, 23, 1211--1216.
  <https://doi.org/10.1016/0009-2509(68)89029-3>
* [UFUA96] Ulshöfer, V. S., Flöck, O. R., Uher, G., and Andreae, M. O.
  (1996). Photochemical production and air-sea exchange of sulfide in the
  eastern Mediterranean Sea. *Mar. Chem.*, 53, 25--39.
  <https://doi.org/10.1016/0304-4203(96)00010-2>

## Diffusivity in soil air and water

* [CH78] Clapp, R. B. and Hornberger, G. M. (1978). Empirical equations
  for some soil hydraulic properties. *Water Resources Res.*, 14(4), 601--604.
  <https://doi.org/10.1029/WR014i004p00601>
* [M03] Moldrup, P., Olesen, T., Komatsu, T, Yoshikawa, S, Schjønning, P, and
  Rolston, D. E. (2003). Modeling diffusion and reaction in soils: X. A
  unifying model for solute and gas diffusivity in unsaturated soil. *Soil
  Sci.*, 168(5), 321--337. <https://doi.org/10.1097/01.ss.0000070907.55992.3c>
"""
module Mass

include("../docstring_style.jl")

using Canopy.Constants: R, atm, T_0

export diffus_air,
    diffus_water,
    diffus_soil_air,
    diffus_soil_water,
    diffus_soil

"""Gas diffusivity in air under STP condition [m^2 s^-1]."""
const diffus_air_stp = Dict{String, Float64}(
    "h2o" => 2.178e-5,
    "co2" => 1.381e-5,
    "ch4" => 1.952e-5,
    "co" => 1.807e-5,
    "so2" => 1.089e-5,
    "o3" => 1.444e-5,
    "nh3" => 1.978e-5,
    "n2o" => 1.436e-5,
    "no" => 1.802e-5,
    "no2" => 1.361e-5,
    "n2" => 1.788e-5,
    "o2" => 1.820e-5,
    "cos" => 1.381e-5 / 1.21)

"""
Pre-exponential factor and activation energy for the calculation of gas
diffusivity in water [m^2 s^-1] under STP condition.
"""
const diffus_water_stp = Dict{String, Tuple{Float64, Float64}}(
    "he" => (818e-9, 11.70e3),
    "ne" => (1608e-9, 14.84e3),
    "kr" => (6393e-9, 20.20e3),
    "xe" => (9007e-9, 21.61e3),
    "rn" => (15877e-9, 23.26e3),
    "h2" => (3338e-9, 16.06e3),
    "ch4" => (3047e-9, 18.36e3),
    "co2" => (5019e-9, 19.51e3),
    "cos" => (4.735872481253359e-6, 19336.20405260121),
    "co" => (0.407e-4, 24518.24),
    "no" => (39.8e-4, 34978.24))

# note: this shall be moved to a soil module once it is started
"""Clapp--Hornberger shape parameters for tortuosity effect."""
const soil_shape_params = Dict{String, Float64}(
    "sand" => 4.05,
    "loamy sand" => 4.38,
    "sandy loam" => 4.9,
    "silt loam" => 5.30,
    "loam" => 5.39,
    "sandy clay loam" => 7.12,
    "silty clay loam" => 7.75,
    "clay loam" => 8.52,
    "sandy clay" => 10.4,
    "silty clay" => 10.4,
    "clay" => 11.4)

"""
Calculate the diffusivity of a gas in air [m^2 s^-1] from temperature [K] and
pressure [Pa].

Support gas species including H2O, CO2, CH4, CO, SO2, O3, NH3, N2O, NO, NO2,
N2, O2, and COS. Note that the `species` argument must be in all lower case.

# See also

[`diffus_soil_air`](@ref)
"""
function diffus_air(species::String, temp, pressure=atm)
    if !haskey(diffus_air_stp, species)
        throw(ArgumentError("species \"$species\" not supported"))
    else
        return diffus_air_stp[species] * (atm / pressure) * (temp / T_0)^1.81
    end
end

"""
Calculate the diffusivity of a gas in water [m^2 s^-1] from temperature [K].

Support gas species including He, Ne, Kr, Xe, Rn, H2, CH4, CO2, CO, NO, and
COS. Note that the `species` argument must be in all lower case.

# See also

[`diffus_soil_water`](@ref)
"""
function diffus_water(species::String, temp)
    if !haskey(diffus_water_stp, species)
        throw(ArgumentError("species \"$species\" not supported"))
    else
        preexp, E_act = diffus_water_stp[species]
        return preexp * exp(-E_act / (R * temp))
    end
end

"""
Calculate the diffusivity of a gas in soil air [m^2 s^-1]. Support gas species
including H2O, CO2, CH4, CO, SO2, O3, NH3, N2O, NO, NO2, N2, O2, and COS.

# Arguments

* `species`: Chemical name of the gas species to be evaluated. In lower case.
* `texture`: Soil texture name in lower case.
* `temp`: Temperature [K].
* `theta_sat`: Total porosity of the soil [m^3 m^-3].
* `theta_w`: Water-filled porosity of the soil [m^3 m^-3].
* `pressure` (*optional*): Ambient pressure [Pa]. Default is the standard
  atmospheric pressure.

# See also

* [`diffus_air`](@ref)
* [`diffus_soil_water`](@ref)
* [`diffus_soil`](@ref)
"""
function diffus_soil_air(species::String, texture::String, temp, theta_sat,
                         theta_w, pressure=atm)
    if !haskey(soil_shape_params, texture)
        throw(ArgumentError("soil texture \"$texture\" not supported"))
    else
        theta_a = theta_sat - theta_w
        tau_a = theta_a * (theta_a / theta_sat) ^
            (3. / soil_shape_params[texture])  # tortuosity in soil air
        return diffus_air(species, temp, pressure) * theta_a * tau_a
    end
end

"""
Calculate the diffusivity of a gas in soil water [m^2 s^-1]. Support gas
species including He, Ne, Kr, Xe, Rn, H2, CH4, CO2, CO, NO, and COS.

# Arguments

* `species`: Chemical name of the gas species to be evaluated. In lower case.
* `texture`: Soil texture name in lower case.
* `temp`: Temperature [K].
* `theta_sat`: Total porosity of the soil [m^3 m^-3].
* `theta_w`: Water-filled porosity of the soil [m^3 m^-3].

# See also

* [`diffus_water`](@ref)
* [`diffus_soil_air`](@ref)
* [`diffus_soil`](@ref)
"""
function diffus_soil_water(species::String, texture::String, temp, theta_sat,
                           theta_w)
    if !haskey(soil_shape_params, texture)
        throw(ArgumentError("soil texture \"$texture\" not supported"))
    else
        tau_w = theta_w * (theta_w / theta_sat) ^
            (soil_shape_params[texture] / 3. - 1.)  # tortuosity in soil water
        return diffus_water(species, temp) * theta_w * tau_w
    end
end

"""
Calculate the total diffusivity of a gas in soil gaseous and aqueous phases
[m^2 s^-1]. Supporting gas species including CO2, CO, NO, CH4, and COS.

# Arguments

* `species`: Chemical name of the gas species to be evaluated. In lower case.
* `texture`: Soil texture name in lower case.
* `temp`: Temperature [K].
* `theta_sat`: Total porosity of the soil [m^3 m^-3].
* `theta_w`: Water-filled porosity of the soil [m^3 m^-3].
* `pressure` (*optional*): Ambient pressure [Pa]. Default is the standard
  atmospheric pressure.

# See also

* [`diffus_soil_air`](@ref)
* [`diffus_soil_water`](@ref)
"""
function diffus_soil(species::String, texture::String, temp, theta_sat,
                     theta_w, pressure=atm)
    if !haskey(soil_shape_params, texture)
        throw(ArgumentError("soil texture \"$texture\" not supported"))
    end

    if !haskey(diffus_air_stp, species) || !haskey(diffus_water_stp, species)
        throw(ArgumentError("species \"$species\" not supported"))
    end

    return diffus_soil_water(species, texture, temp, theta_sat, theta_w) +
        diffus_soil_air(species, texture, temp, theta_sat, theta_w, pressure)
end

end  # module


text:
# Serve WGLMakie via web

import JSServe
JSServe.configure_server!(listen_url="0.0.0.0", external_url="http://192.168.0.101:8082", listen_port=8082) # does not work
# JSServe.configure_server!(listen_port=8082) # works
using WGLMakie
WGLMakie.activate!()

fig = Figure()

ax = Axis(fig[1, 1])

lsgrid = labelslidergrid!(
    fig,
    ["Voltage", "Current", "Resistance"],
    [0:0.1:10, 0:0.1:20, 0:0.1:30];
    formats = [x -> "$(round(x, digits = 1))$s" for s in ["V", "A", "Ω"]],
    width = 350,
    tellheight = false)

fig[1, 2] = lsgrid.layout

sliderobservables = [s.value for s in lsgrid.sliders]
bars = lift(sliderobservables...) do slvalues...
    [slvalues...]
end

barplot!(ax, bars, color = [:yellow, :orange, :red])
ylims!(ax, 0, 30)

set_close_to!(lsgrid.sliders[1], 5.3)
set_close_to!(lsgrid.sliders[2], 10.2)
set_close_to!(lsgrid.sliders[3], 15.9)

display(fig)

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
export marginalrule

@marginalrule MvNormalMeanPrecision(:out_μ_Λ) (m_out::MultivariateNormalDistributionsFamily, m_μ::PointMass, m_Λ::PointMass) = begin
    return (out = prod(ProdAnalytical(), MvNormalMeanPrecision(mean(m_μ), mean(m_Λ)), m_out), μ = m_μ, Λ = m_Λ)
end

@marginalrule MvNormalMeanPrecision(:out_μ_Λ) (m_out::PointMass, m_μ::MultivariateNormalDistributionsFamily, m_Λ::PointMass) = begin
    return (out = m_out, μ = prod(ProdAnalytical(), m_μ, MvNormalMeanPrecision(mean(m_out), mean(m_Λ))), Λ = m_Λ)
end

@marginalrule MvNormalMeanPrecision(:out_μ_Λ) (m_out::MultivariateNormalDistributionsFamily, m_μ::MultivariateNormalDistributionsFamily, m_Λ::PointMass) = begin
    xi_y, W_y  = weightedmean_precision(m_out)
    xi_m, W_m  = weightedmean_precision(m_μ)

    W_bar = mean(m_Λ)

    T = promote_type(eltype(W_bar), eltype(W_y), eltype(W_m))
    d = length(xi_y)
    Λ = Matrix{T}(undef, (2*d, 2*d))
    @inbounds for k2 in 1:d
        @inbounds for k1 in 1:d
            tmp1 = W_bar[k1,k2]
            tmp2 = -tmp1
            k1d = k1+d
            k2d = k2+d
            Λ[k1,k2] = tmp1 + W_y[k1,k2]
            Λ[k1d,k2] = tmp2
            Λ[k1,k2d] = tmp2
            Λ[k1d,k2d] = tmp1 + W_m[k1,k2]
        end
    end

    ξ  = [ xi_y; xi_m ]

    # naive:
    # Λ  = [ W_y + W_bar -W_bar; -W_bar W_m + W_bar ]
    # ξ  = [ xi_y; xi_m ]

    return (out_μ = MvNormalWeightedMeanPrecision(ξ, Λ), Λ = m_Λ)
end

@marginalrule MvNormalMeanPrecision(:out_μ) (m_out::MultivariateNormalDistributionsFamily, m_μ::MultivariateNormalDistributionsFamily, q_Λ::Any) = begin
    xi_y, W_y  = weightedmean_precision(m_out)
    xi_m, W_m  = weightedmean_precision(m_μ)

    W_bar = mean(q_Λ)

    T = promote_type(eltype(W_bar), eltype(W_y), eltype(W_m))
    d = length(xi_y)
    Λ = Matrix{T}(undef, (2*d, 2*d))
    @inbounds for k2 in 1:d
        @inbounds for k1 in 1:d
            tmp1 = W_bar[k1,k2]
            tmp2 = -tmp1
            k1d = k1+d
            k2d = k2+d
            Λ[k1,k2] = tmp1 + W_y[k1,k2]
            Λ[k1d,k2] = tmp2
            Λ[k1,k2d] = tmp2
            Λ[k1d,k2d] = tmp1 + W_m[k1,k2]
        end
    end

    ξ  = [ xi_y; xi_m ]

    # naive:
    # Λ  = [ W_y + W_bar -W_bar; -W_bar W_m + W_bar ]
    # ξ  = [ xi_y; xi_m ]

    return MvNormalWeightedMeanPrecision(ξ, Λ)
end

@marginalrule MvNormalMeanPrecision(:out_μ) (m_out::PointMass, m_μ::MultivariateNormalDistributionsFamily, q_Λ::Any) = begin
    return (out = m_out, μ = prod(ProdAnalytical(), MvNormalMeanPrecision(mean(m_out), mean(q_Λ)), m_μ))
end

@marginalrule MvNormalMeanPrecision(:out_μ) (m_out::MultivariateNormalDistributionsFamily, m_μ::PointMass, q_Λ::Any) = begin
    return (out = prod(ProdAnalytical(), MvNormalMeanPrecision(mean(m_μ), mean(q_Λ)), m_out), μ = m_μ)
end

## ProdFinal / BIFM related

@marginalrule MvNormalMeanPrecision(:out_μ_Λ) (m_out::ProdFinal, m_μ::PointMass, m_Λ::PointMass) = begin
    return (out = getdist(m_out), μ = m_μ, Λ = m_Λ)
end

@marginalrule MvNormalMeanPrecision(:out_μ_Λ) (m_out::PointMass, m_μ::ProdFinal, m_Λ::PointMass) = begin
    return (out = m_out, μ = getdist(m_μ), Λ = m_Λ)
end


text:
using Test, PropDicts

@testset "Get Set" begin
    d = PropDict(a=1, b=2)
    @testset "get" begin
        @test d.a == 1
        @test d.b == 2
        # unknown key
        @test_throws KeyError d.k
    end

    @testset "set" begin
        d.a = 7
        @test d.a == 7

        # unknown key no wuzzas
        d.j = "11"
        @test d.j == "11"
    end

    @testset "typed" begin
        dt = PropDict{Int}(a=1)
        @test @inferred getproperty(dt, :a) == 1
        @test_throws ErrorException @inferred getproperty(d, :a)
    end
end


text:
# Ref:
# https://en.wikipedia.org/wiki/Inside%E2%80%93outside%E2%80%93beginning_(tagging)
# https://chameleonmetadata.com/Education/NLP-3/ref_nlp_encoding_schemes_list.php

abstract type tag_scheme end

struct BIO1 <: tag_scheme end # BIO
struct BIO2 <: tag_scheme end
struct BIOES <: tag_scheme end

const available_schemes = ["BIO1", "BIO2", "BIOES"]

"""
    tag_scheme!(tags, current_scheme::String, new_scheme::String)

Convert `tags` from `current_scheme` to `new_scheme`.

List of tagging schemes currently supported-
 * BIO1 (BIO)
 * BIO2
 * BIOES

# Example
```julia-repl
julia> tags = ["I-LOC", "O", "I-PER", "B-MISC", "I-MISC", "B-PER", "I-PER", "I-PER"]

julia> tag_scheme!(tags, "BIO1", "BIOES")

julia> tags
8-element Array{String,1}:
 "S-LOC"
 "O"
 "S-PER"
 "B-MISC"
 "E-MISC"
 "B-PER"
 "I-PER"
 "E-PER"
```
"""
function tag_scheme!(tags, current_scheme::String, new_scheme::String)
    current_scheme = uppercase(current_scheme)
    new_scheme = uppercase(new_scheme)
    (length(tags) == 0 || current_scheme == new_scheme) && return

    if new_scheme ∉ available_schemes || current_scheme ∉ available_schemes
        error("Invalid tagging scheme")
    end

    current_scheme = eval(Symbol(current_scheme))()
    new_scheme = eval(Symbol(new_scheme))()

    tag_scheme!(tags, current_scheme, new_scheme)
end

function tag_scheme!(tags, current_scheme::BIO1, new_scheme::BIO2)
    for i in eachindex(tags)
        if tags[i] == 'O' || tags[i][1] == "O"
            tags[i] = "O"
            continue
        end
        (tags[i][1] == 'O' || tags[i][1] == 'B') && continue

        if tags[i][1] == 'I'
            if i == 1
                tags[i] = 'B' * tags[i][2:end]
            elseif tags[i - 1] == "O" || tags[i - 1][2:end] != tags[i][2:end]
                tags[i] = 'B' * tags[i][2:end]
            else
                continue
            end
        else
            error("Invalid tags")
        end
    end
end

function tag_scheme!(tags, current_scheme::BIO2, new_scheme::BIO1)
    for i in eachindex(tags)
        if tags[i] == 'O' || tags[i][1] == "O"
            tags[i] = "O"
            continue
        end
        (tags[i][1] == 'O' || tags[i][1] == 'I') && continue

        if tags[i][1] == 'B'
            if i == length(tags)
                tags[i] = 'I' * tags[i][2:end]
            elseif tags[i + 1] == "O" || tags[i + 1][2:end] != tags[i][2:end]
                tags[i] = 'I' * tags[i][2:end]
            else
                continue
            end
        else
            error("Invalid tags")
        end
    end
end

function tag_scheme!(tags, current_scheme::BIO2, new_scheme::BIOES)
    for i in eachindex(tags)
        if tags[i] == 'O' || tags[i][1] == 'O'
            tags[i] = "O"
            continue
        end

        if tags[i][1] == 'I' && (i == length(tags) ||
                                 tags[i+1][2:end] != tags[i][2:end])
            tags[i] = 'E' * tags[i][2:end]
        elseif tags[i][1] == 'B' && (i == length(tags) ||
                                 tags[i+1][2:end] != tags[i][2:end])
            tags[i] = 'S' * tags[i][2:end]
        else
            (tags[i][1] == 'I' || tags[i][1] == 'B') && continue
            error("Invalid tags")
        end
    end
end

function tag_scheme!(tags, current_scheme::BIOES, new_scheme::BIO2)
    for i in eachindex(tags)
        if tags[i] == 'O' || tags[i][1] == 'O'
            tags[i] = "O"
            continue
        end
        (tags[i][1] == 'B' || tags[i][1] == 'I') && continue

        if tags[i][1] == 'E'
            tags[i] = 'I' * tags[i][2:end]
        elseif tags[i][1] == 'S'
            tags[i] = 'B' * tags[i][2:end]
        else
            error("Invalid tags")
        end
    end
end

function tag_scheme!(tags, current_scheme::BIO1, new_scheme::BIOES)
    tag_scheme!(tags, BIO1(), BIO2())
    tag_scheme!(tags, BIO2(), BIOES())
end

function tag_scheme!(tags, current_scheme::BIOES, new_scheme::BIO1)
    tag_scheme!(tags, BIOES(), BIO2())
    tag_scheme!(tags, BIO2(), BIO1())
end